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Summary 


The electrical analogue method of investigating heat transfer processes was used 


to determine the thermal characteristics of proposed wall and roof structures for a 40,000 ton 


raw sugar silo. 
consequent transport of moisture. 


Thermal insulation is necessary to prevent internal temperature gradients and 
Simplification of the problem was achieved by separating 


the external temperature factors into a daily variation of air temperature, a seasonal variation 


and solar radiation. 
one of the analogue networks. 


The details of the method are illustrated by calculating the constants of 
Results for seven structures are tabulated. 


Effective thermal conductivities of 19 packed raw sugars were measured to provide data not 


previously available. 
the bulk density of the packed sugar. 


Résumé 


The conductivities are shown to be related to the moisture content and 


Pour déterminer les caractéristiques, thermiques des murs et de la toiture d'un silo a 


sucre brut d'une capacité de 40.000 tonnes, les auteurs utilisent la méthode d’étude des processus 


de transfer de chaleur dite 


“méthode d’analogie électrique’’. 


Une isolation thermique permet d'éviter des gradients de température interne et par conséquent 


le transport dhumidité. 


Le probléme se simplifie en séparant les facteurs externes de la tem- 


pérature en une variation journaliére de la température de l'air, une variation saisonniére et une 


irradiation solaire. 


Les auteurs illustrent les détails de la méthode par le calcul des constantes d'un des réseaux 


d'analogue. 


Ils donnent leurs résultats pour le calcul de sept constructions. 


Ils mesurent les 


conductivités thermiques effectives de dix-neuf empilements de sucre brut afin d’obtenir des 


données manquant jusqu’é présent. 
occluse et & la densité brute du sucre empilé. 


INTRODUCTION 


When raw sugar is stored in bulk there is a 
possibility of loss of sugar on account of the 
activity of micro-organisms (osmophilic yeasts 
and moulds) which increases with increasing 
moisture content. Raw sugar contains a small 
quantity of water, of the order of 05%. The 
percentage is variable but is usually less than 
the value at which appreciable destruction of 
sugar occurs. If temperature gradients were to 
arise within the sugar they would cause thermal 


Ils montrent que les conductivités sont liées & humidité 


transport of moisture towards the regions at the 
lower temperatures and the moisture content 
might be increased locally above the threshold 
Furthermore, the regions at the higher 
temperatures would be hardened by loss of 
moisture and difficulties would be experienced 
when the sugar was reclaimed from the silo. 


value. 


A silo designed to hold 40,000 tons of raw sugar 
was recently built at the Plaistow Wharf (London) 
Refinery of Tate and Lyle Ltd. [1]. The silo has 
a retaining wall in the form of a vertical steel 
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cylinder of diameter 160 ft. and height 90 ft. and 
Full 
capacity is achieved with a depth of sugar of 
80 ft. 
circular cylinder of fixed volume resting on the 


it is spanned by a slightly domed roof. 


The total exposed surface area of a solid 


ground with its axis vertical is a minimum when 
its height is equal to its radius. This condition is 
fulfilled by the cylinder of sugar at full capacity 
80 ft.) but in spite of this the 
and the 
mass of sugar within 8 ft. of the surface is about 
6,000 tons. 


sugar in this boundary 


(height radius 


surface area amounts to 60,000 sq ft. 


The condition during storage of the 
laver is a matter of 
considerable importance. 

Temperature gradients near to the surface of 
the sugar due to changes of ambient temperature 
must therefore be reduced by the provision of 
adequate thermal insulation of the walls and roof 
of the 
proposed 


silo. Certaim insulating systems were 


and their thermal properties were 


investigated by means of the electrical analogue 
method [2, 3] 
ANALOGUI 


THEORY Oo} ELECTRICA! 


The 


system of 


rut 


distribution of temperature throughout a 
heat 


media in contact with each other is determined 


uniform, isotropic conducting 


at all pomts by Fourier’s equation : 


(1) 


Similarly the distribution of potential in an 
electrically conducting system is determined by 
an analogous equation : 


vv=-Rc.? 


. 


in which R is the resistivity and C is the capaci- 
tance per unit volume. 

p, c and k in the thermal case, and R and C in 
the electrical case are constant within each region 
of the system but vary from one region to 
another. 

The electrical analogue method is based on the 
similarity of equations (1) and (2). 

Analogous quantities are : 


Temperature : 
Heat : 

Flow of heat 
Time : 


Voltage 

Electric charge 
Electric current 
Time 


Thermal resistance : Electrical resistance 


Thermal capacity : 


Electrical capacitance 


Between any points of the electrical analogue 
the values of electrical resistance and capacitance 
are made proportional to the values of their 
thermal counterparts between the corresponding 
points of the thermal system. 
the 
analogue a practical proposition it is necessary 


In order to make construction of an 
to use lumped resistors and capacitors in place 
of the theoretically ideal distributed impedances. 
The analogue now consists of a two or three- 
dimensional network of resistors with capacitors 
connected between earth and the junction points 
of the resistors. Lawson and McGuire [4] have 
discussed the magnitude of the errors involved 
by this substitution. Actual physical dimensions 
of the analogue units are, of course, of no signifi- 
cance. Corresponding points of the thermal and 
electrical systems are determined by the resist- 
An 
important advantage of the analogue method is 
Whereas the 
time constant of a section of the thermal system 


ances and capacitances between the points. 
the reduction of the time scale. 


might be of the order of days, that of the corres- 
ponding section of the analogue may be reduced 
if necessary to the order of microseconds. Cathode 
ray tube display of the voltage variations ts a 
convenient means of indication. Temperature 
variations at the boundary of the thermal system 
the 
corresponding points of the analogue. These have 
but 


respect to the temperature variations by the 


are simulated by voltage variations at 


a similar wave form are accelerated with 
appropriate time scale ratio. Generation of heat 
at a point in the thermal system and flow of heat 
into it from outside are simulated by flow of 
current, of the the 


corresponding points of the analogue. 


appropriate value, into 


APPLICATION TO THE Sito PROBLEM 


The radii of curvature of the walls and roof of 
the silo are very large compared to the depth of 
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penetration of the temperature variations into 
the sugar, so that for the purposes of this problem 
the walls and roof can be considered to be plane. 
Heat flow in the sugar and in the insulating 
layers is normal to the surface and is one dimen- 
sional. Thus the electrical analogue becomes a 
single line of resistors with capacitors to earth 
(Fig. 1). 


aseestos asetsTos 


wremrace 


Fic. 1. 


Ist design, for daily variations of air temperature. 
ances in thousands of ohms. ¢ 


Z 


To simplify the problem the three sources of 


temperature variation were considered separately. 


The sources are : 


Daily variations of temperature of the air 
outside the silo, 
Seasonal variations of temperature of the air 
outside the silo, 
Solar radiation. 


Since rate of heat transfer is proportional to 
the first power of temperature difference the 
measured 


effects of these three factors can be 


separately and added together. This is only 
approximately true for transfer across an air 
gap. Each of the six silo wall structures that 
were investigated included an air gap, but, as 
the temperature differences across the gaps are 
small, the approximation was considered to be 
sufficiently accurate. 

Records show that 
between maximum and minimum temperatures 
within a 24-hour period is not usually greater 
than 14°C and that this is also the difference 
between average July and January temperatures. 
The daily and seasonal variations can be treated 


in London the difference 


measurement by electrical analogue 


with sufficient accuracy as sine waves of amplitude 
7°C. The 
penetrating power than the fundamentals on 


harmonics are small and have less 
account of their higher frequencies. This approxi- 
mation is justified by the oscillograms shown in 
Fig: 5. Both variations were represented at the 
input to the analogue by a 50 evcle sec. voltage 


from the A.C. mains, so that two analogues, each 
JOSCKLOK OME 
ad ae 


woow> 


ee eo 


eo 
Wy AA Oe 
| cts 

a = = —_ —_— — lz 
K ss pose lo 0425 00425 fosrs boas 
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Electrical analogue of the bottom of the silo wall, 


Resist- 
apacitances in microfarads. 


characteristic impedance of the sugar sections. 


with its appropriate time constants, were required 
for each insulating system. The analogue of solar 
radiation is a current flowing into the input end 
of the electrical line having a wave form corres- 
ponding to the daily variation of the flux of solar 
radiation. The form of this variation changes 
throughout the year, and the most severe con- 
ditions should be investigated if possible. On the 
roof most severe conditions occur on a cloudless 
day in midsummer. Consideration of the variation 
of severity of conditions affecting the wall is very 
complex. It was considered sufficient to investi- 
gate midsummer conditions for the south facing 
part of the wall. 

It is shown in Appendix 2 that, in midsummer 
at latitude ¥, the intensities of radiation on a 
horizontal surface (roof) and on the south facing 
part of the silo wall are respectively : 

Ne = No (0-92 cos w sin wt + 0-40 sin p) 
and 


Nw Ng (0-92 sin % cos wt + 0-40 cos Wp) 


These expressions are applicable when they are 
positive. Negative values of Np, mean that the 
intensity is zero as the sun is below the horizon. 
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FIBREGLASS 


< — 
ASBESTOS . 
; 


4 


OUTSIDE 











Section through the bottom of the silo wall, 
design. 





iw 
wt 
Fic. 3. 


Incidence of solar radiation. 


When Ny is negative the sun is behind the plane 
of a south facing wall. 


London is at latitude 514° so that 


Ne 
and 
Ny Ny (0-72 cos wil 


No (0-57 cos wil 0-31) 


0-25) 
These wave forms are shown in Fig. 3 (a and b). 


Interrupted cosine current waves required to 


and Ny 
produced by means of the circuit shown in Fig. 4, 


satisfy the expressions for N, were 
in which the valve is biassed so that the correct 
proportion of the full wave is suppressed. Fig. 5(a) 
shows an oscillogram of the input current wave 
form. 

The effects of the three sources of t mperature 
could have been obtained by direct 
but the 


involved the lengthy solution of sets of simul- 


Variation 


calculation calculation would have 


taneous equations derived from the boundary 
conditions at each of the interfaces between the 


layers. In addition it would have been necessary 


MULLER 


of solar 


effect 


radiation into its Fourier components and solve 


to resolve the non-sinusoidal 
the equations for each component separately. 
Furthermore each of the systems investigated 
would have required a complete calculation, 
the and 
measured in a short time. 


whereas analogues were assembled 


CALCULATION OF THE 


CONSTANTS OF AN ANALOGUE 


The procedure for detailed calculation of an 
analogue will now be illustrated by reference to 
the “ daily ~ 
proposed for the silo wall. The layer of steel which 
bears the lateral thrust of the sugar tapers from 


analogue of one of the systems 
‘ b 


' 
0-72 
\ wr } iv 


\ O25 
= | ’ 








(a) on the roof; (b) on the south wall. 





— 


TO ANALOCUE 

















Fic. 4. Circuit for interrupted cosine waves of current. 


bottom to top. The analogue given here repre- 
sents the bottom of the wall where the thickness 
of the shell is greatest. A full specification of the 
bottom of the wall is given in Table 1 below. 
For the purpose of heat transmission the wall 
consists of the five layers indicated in Fig. 2, 
which is not to scale. Their relevant properties, 
together with those of an average raw sugar are : 
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Table 1. 





Air 
asbestos 
interface 


MATERIAL 


THICKNESS (metres) 0 


Density (kg/metre*) 
Mass /metre* 
(kg /metre*) 


Se. Hr. 
(joules /kg °C) 


THERMAL CONDUCTIVITY 
(watts ‘metre “C) 

THERMAL RESISTANCE 0-053 
per metre* (°C /watt) 

THermMaL CAPACITY 

per metre? (joules /°C) (=c xX m) 


Time CONSTANT 
(seconds) 


Asbestos 


1762 80 


10350 


Fibre- 
glass 


Air gap 


0-00635 0-0762 0-381 0-0203 


7800 910 


| 910/ 
| metre depth 


0-033 


0-023 

metre depth 

| 1-14 x 10% 

| metre depth 

| 71x 108 
metre” 





Note : 1 day = 86,400 seconds. 


The thermal resistances of the air-solid interface 
and the air gap were taken from The Computation 
of Heat Requirements for Buildings, published by 
the Institution of Heating Ventilating 
Engineers, 1950. 

In row 9 of Table 1 the “time constant” is 
the product of the thermal resistance and thermal 
capacity of a square metre of the particular layer 
It has the dimensions of time only 


and 


considered. 

and is independent of the cross sectional area. 
Let us select any one of the layers that has both 

We 


resistance and capacity — say the fibreglass. 


must satisfy the condition :; 
Time constant of the fibreglass section of the analogue 
Time scale of the electrical signal 
Time constant of the fibreglass 


Time scale of the temperature variation 
9,300 seconds 


: —— = 0:108 
86,400 seconds 


The period of the electrical signal from the mains 
is 0-02 seconds. Hence the time constant of the 
fibreglass section 


0-108 « 0-02 secs. = 2-16 x 10° secs. 


We may now arbitrarily set the capacitance of 
the fibreglass section to some convenient value 


say 0-003 microfarads — and enter this value* in 


row 8 of Table 2. 


Then since RC = 2-16 10° sees. 


9.18 3 
R 2-16 «x 10 — 7-2 
38 x 10° 


10° ohms 


This value’ is entered in row 2 of Table 2. The 
other sections of the analogue are given resistances 
and capacitances which are proportional to their 
thermal resistances and capacities, the values 
being based on 0-003 » farads and 7-2 x 10° ohms 
for the fibreglass section. 

EXPERIMENTAL PROCEDURE 
(a) Daily and Seasonal Temperature Variations 


Each of the analogues was assembled as shown 
in Fig. 1. The amplifier of the cathode ray oscillo- 
graph was connected between earth and a probe 
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Table 2. 





Air 
MATERIAI asbestos Asbestos Fibre- Air gap Sugar 
interface glass 


ELECTRICAL RESISTANCE 16400 T20000F 84000 19200 com depth 


(ohms) 


ELECTRICAL CAPACITANCE ; 0-0078 0-00380* 0-00035 0-053 0-0085, em depth 


(nu farads) 








Ideally the network should be terminated as shown by an impedance equal to the characteristic impedance of the 
sugar sections, but this termination was not actually used as the amplitude of the signal was reduced almost to zero at 
the fifth sugar section so that spurious reflections were negligible. 

In the “ seasonal’ analogue each of the five sugar sections represented a 100 em depth of sugar 


Table 8. 





Insulation factor 


Structure Materials Thickness of 
layer (cm) Daily Seasonal 


Top of wall (ist design) Asbestos 0-635 14 
Fibreglass 7 é2 
Air gap 38-1 
Steel 0-38 


Bottom of wall (Ist design) Asbestos 0-635 
Fibreglass 7 62 
Air gap 


Steel 


Top of wall (2nd design) Asbestos 
Air gap 
Steel 


Bottom of wall (2nd design) Asbestos 
Air gap 
Steel 


Roof 1 Foam concrete 


Roof 2 Ruberoid 0-635 
Steel 317 


Top of wall of Swedish silo Steel OS 
Air gap 15-2 
Wood 3-81 


Bottom of wall of Swedish silo Steel 2-03 
Air gap 15-2 
Wood 3-81 





Design 2 differs from design 1 in that it has no fibreglass. An air space 6 metres thick was assumed to exist below 


the roof 
The Swedish silo above represents an existing structure at Gothenburg and is included here for comparison. 
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(ce) 


Fic. 5. Temperature wave forms due to solar radiation. 

(a) at the outer surface of the asbestos (B); (b) at the 

outer surface of the sugar (EF); (c) 125m inside the 
sugar (F'). Refer to Fig. 1 for points B, E and F. 








Thermal insulation of a raw sugar silo : 


which could be applied to each of the points A, 
»C. 


amplifier gain to give a predetermined amplitude 


Having set the input voltage and the 


on the oscillograph screen with the probe at the 
point A, the probe was applied to the points B, 
i w« 
The 


given by an “ 


and the amplitudes were observed. 
the 


insulation factor ” 


effectiveness of thermal insulation is 
which is defined 
as the ratio of the amplitudes at the 
A and E of the system. At 


any of the internal points B, C, 


points 


D...the phase of the signal is 
changed with respect to the input 
This phase change can be 
the 
sweep of the oscillograph is locked 
the 50 


signal. 


readily observed if horizontal 


directly with ye le see 


supply (Fig. 5). 


(b) Solar Radiation 


Having decided, arbitrarily, the re- 


lation between voltage and tem- 


the 


represent a given flux of radiation 


perature current required to 
can be calculated by noting that the 


dimensionless quantities : 


Temperature x Capacity /unit area 


Time scale Heat flux unit area 
Voltage Capacitance 
Time scale Current 


are equal, 


Heat transferred from the surface 


to the atmosphere by conduction 


and convection must be allowed for 
by connecting the resistance R of 
Fig. 1, 
film resistance, from the point B to 
earth. When measuring the effect of 


the analogue of the surface 


radiation on the south wall R was 
increased by a factor 5/3 to take 
account of the fact that the surface 
resistance of a 5/3 
times the average for the wall as a 


Fic. 6. 
wall is 


south 


whole. For the roof the surface film 
resistance was taken as half of that 


WALL OF SWEDISH SILO 
(TOP AND BOTTOM) 


Temperature ranges (range 2 

Column 1. 

Column 2. 

Column 3. 

The height of the left hand side of a shaded area represents the range of 

temperature at the outer surface of the wall or roof; the height at the 
right hand side represents the range at the surface of the sugar. 


measurement by electrical analogue 


for a south wall. A current, adjusted to the correct 
wave form by the setting of R and V in the 
circuit of Fig. 4 is fed in at B. 

As the input is not now sinusoidal the wave 
form of the voltage-time curve changes from 
point to point of the network. Oscillograms (which 
are relatively not to scale) of some of the wave 


forms are shown in Fig. 5. It can be seen that 


DAILY 
VARIATION 


SEASONAL 
VARIATION 


OUTER SURFACE 
INSULATION 


INSULATION 
INSULATION 


< amplitude). 
Due to daily temperature variation. 

Due to seasonal temperature variation. 
Due to solar radiation in midsummer. 
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the wave form becomes more nearly sinusoidal 


as it penetrates more deeply into the system. 


RESULTS 

Details of the proposed wall and roof structures 

which were investigated are Table 8, in 

the the 
Insulation factors are given in columns 


NUMERICAL 


given in 


which first mentioned layers are on 
outside, 
4 and 5. 

The insulation factors are illustrated graphic- 
The heights of the left and right 


indicate the ranges 


ally in Fig. 6. 
hand side of a shaded area 
of temperature at the outer surface of the wall 
and the boundary of the sugar respectively, An 
absorption coeflicient of 0-5 has been assumed for 
the effective 


radiation. The solar radiation diagrams re present 


the caleulation of flux of solar 
midsummer conditions. 


Two points of gene ral interest can be noted. 


(1) In all cases imsulation is much less effective 
the effect of the 


variation of temperature than in suppressing 


in suppressing NO asonal 
the daily variation. This is due to the effect 
of thermal capacity. 


The two silo walls (Ist and 2nd designs) have 
greater insulation factors at the bottom than 
at the top. This difference 
between the thermal capacities per unit area 
of the steel shell at the top, where the thick- 
ness is 0-38 cm, and at the bottom, 
2-03 cm. The 
appreciable thermal resistance, even at its 
thickest part. 


is due to the 


wihiere the 


steel has 


thickness is 


APPENDIX I! 
As data.on the effective thermal conductivities of raw 
sugars when packed under normal compressions were not 


available, measurements were made on nineteen raw 


sugars. Percentage moisture and bulk density were also 
observed. 


Table 4. 
For the purpose of designing the analogues, values 


Results are shown in 


k 0-161 watts /metre °C 
and p 0-910 kg /litre 


were taken as representative of an average raw sugar. 
It is reasonable to suppose that the conductivity 


increases with moisture and with bulk density. If we 


MULLER 


Table 4. 





Origin of Moisture 
raw sugar (“o) 


Density 
(kg /litre) 


Conductivity 
(watts , metre 
C.) 
0-190 
0-182 


0-265 








British West 
Indies 


0-903 
0-910 
0-857 





0-940 
0-860 
0-936 
0-879 


0-188 
0-157 
174 
O-155 


San 
Domingo 





0-930 


- 
OuUnT 


154 
O-171 
O-l42 


Australia 
OSES 





S79 
928 
‘O4 


O14 
0-169 
156 


Mauritius 





S74 0-165 
v4Y 0-198 
sou 0-173 


England 
(beet sugar) 





910 
SRT 0-152 


825 Os! 


0-160 





O016 00-1642 





ussume a linear dependence of conductivity on moisture 


and density the regression equation is : 


k = 0-0527 0-370 p — 01925 


This equation holds within the limits : 


Ol4< 0-85 and 0-825 < p < 0.970 
Statistical analysis shows that the coefficients of W 


and p are both significant. The standard deviation of 
the measured values of & from the values given by the 
equation is 0-0086, giving a coefficient of variation of 5%. 
If its 
moisture had been the average value W the conductivity, 
been 0-0527 
If we adjust all 
of the conductivities in this way we have a set of con- 


Consider a raw sugar having moitsure W%,. 


according to the equation, would have 


(Ww WW) less than the observed value. 


ductivities corrected for non-average moisture which we 
can plot, as in Fig. 7 (a) against bulk density. 

In Fig. 7 (b) the conductivities, adiusted for non-average 
density, are plotted against moisture. 


In both cases the correlation is good. 
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CONDUCTIVITY 





BB 92 


DENSITY 


—_—e 


(a) 
Fic. 7. 


measurement by electrical analogue 





CONDUCTIVITY — = watts/m °C 








0-4 
MOISTURE 


(hb) 


Thermal conductivities of raw sugar in bulk. 


(a) Conductivity adjusted for non-average moisture plotted against bulk density. (b) Conductivity adjusted for 
nonaverage density plotted against moisture. 


APPENDIN II 
Calculation of the sun's elevation 


Take the rectangular axes Orv, Oy, and Oz fixed in the 
earth and rotating with it (Fig. 8) 


Z 


TO LONDON 


EQUATOR 


Fic. 8 


Ox and Oy are in the equatorial plane. 

Ox passes through the meridian of the locality in 
In midsummer the angle 
234°. 


question, Oz is the polar axis. 
between OS and the Ory plane yp’ 
The angle between OL and the Oxy plane = ¢ 
Let the direction cosines of OS be /,, m,, n, and of 
OL be Ig, ma, Ne- 


514°. 


Then iF cos J’ cos al, cos ob 
my cos d’ sin wf. 0 


ur sin sin ob 


Now the angle of incidence of radiation on the silo roof 


LOS ¢ 
cos } Il, 
cos ¢' cos % cos wt + sin f’ sin & 
0-57 cos wt + O31 


mM, Me + M No 


The south facing wall of the silo is parallel to the earth's 
surface in latitude (¢ - 
Hence the direction cosines ure 


w/2). 


l, = sin ¢ 
ms, = 0 
ng - cos % 


If @ angle of incidence on the wall : 


i 


cos 0 = cos wf’ sin ¢ cos wt — sin ’ cos yp 


= (72 cus wl — 0-25 


APPENDIN III 
Propagation of Temperature Waves into Packed Raw Sugar 
Consider a semi-infinite mass of packed raw sugar on one 
side of a plane boundary. At any instant the temperature 
of the boundary is uniform. The temperature varies with 
time according to : 


T = Ty cos wi (1) 


Heat flow in the sugar is perpendicular to the boundary 
and obeys the equation : 
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The solution of equation (2) which satisfies the boundary ACKNOWLEDGEMENTS 


condition (1) is : : 
— Acknowledgment is made to E. T. Moss for 


T « Tp con wt ! | -2y (wii information regarding the thermal properties of 


| V (2a/w)] the materials considered, to H. C. S. pe WHALLEY, 
Director of Research, and to the Directors of 
Tate & Lyle Ltd. for permission to publish this 


2a 
A 2», / - with an attenuation coefficient 8 paper. 


This equation represents a wave of wavelength 


The relevant constants for raw sugar are NOTATION 
4 temperature 
p = 910 kg/ metre’ 
5 900 feulen fen .y, = = rectangular coordinates 
o« “a “~— K 
. ; & density 
0-161 watts metre “¢ - 
c = specific heat 
Hence « 2» 7 metres? sec thermal conductivity 
time 
~ 
Daily Temperature Variation wv dy? 
, . : thermal diffusivity 
hours 7-26 X 10% sex electric potential 


2r 
> resistance 
So that A Bll ‘ = capacitance 
io = time constant 
0-40 metres , intensity of radiation on the roof 
[: i. } intensity of radiation on the south wall 


\ 2-84 solar constant 


16 metres! 980 watts metre® in London allowing for atmos- 


pheric absorption 
Seasonal Temperature Variations : thermal resistance 
- thermal capacity 
27 ™ ) 
hours’ =< 1-09 « 10°7 sec"! mass per unit area of a layer 
8766 3 
thickness of a layer 
So that A 7-5 metres % moisture 
' , 2 & frequency 
wavelength 
The values of 8 are such that the amplitude of the daily 8 = attenuation coefficient 
temperature wave decreases to 10°, of its boundary value latitude 514° for London 
after a penetration of l4cm; that of the seasonal wave y latitude of the tropics 234 
after a penetration of 2-8 metres. .m, direction cosines 
These figures, together with the wavelengths, were angle of incidence on the roof 
angle of incidence on the south wall 


8 O85 metres 


confirmed by the analogue measurements. 
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Stationary outer cylinder and a rotating inner cylinder 
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(Received 12 February 1955) 


Summary 


The longitudinal mixing produced in a liquid which flows through the space between 


a stationary outer cylinder and a rotating cylindrical core can be interpreted as the result of an 


apparent longitudinal diffusivity. 


For various dimensions and rotor speeds the value of the 


diffusion coefficient has been measured by means of the frequency response analysis. 


Résumé 


une rotation du cylindre intérieur produit un mélange en direction axiale. 
étre interprété comme le résultat d'un coefficient apparent de diffusion longitudinale. 


Dans un liquide, qui se trouve dans espace annulaire entre deux parois cylindriques, 


Ce phénoméne peut 
La valeur 


de ce coefficient a été determiné dans des conditions diverses au moyen de lanalyse de réponse 


en fonction de la fréquence. 


INTRODUCTION 


In a previous article [1] attention has been drawn 
to the use of frequency response analysis as a 
means of determining the distribution of resi- 
dence times. This analysis consists of introducing 
a sinusoidally varying signal in the flow at the 
inlet of the flow system and of measuring the 
amplitude ratio and the phase shift between the 
signals at two or more places along the main flow 
path at various frequencies. 

The result of such an analysis can be compared 
with the calculated behaviour of relatively simple 
models such as a “perfect” mixer, a cascade of 
perfect mixers, or “perfect” piston flow. The 
latter shows no damping of the signal but only 
a phase shift caused by the distance-velocity lag. 
This kind of flow is approached in the turbulent 
flow of a fluid through a straight length of pipe 
and through packed beds, but a still significant 
deviation may occur as a consequence of mixing in 
the direction of flow. If this mixing is of a suffi- 
ciently random nature, it can be accounted for 
mathematically by assigning an apparent co- 
efficient of diffusion, or a longitudinal diffusivity, 
to the fluid stream, Thus, if longitudinal diffusion 
occurs in a flow system, the diffusivity can be 
determined by carrying out a frequency response 


analysis. 


The present paper reports on such measure- 
ments for the flow of a liquid through the annular 
space between a rotating inner cylinder and a 
concentric stationary outer cylinder. This parti- 
cular choice was mainly determined by the fact 
that the 
depends primarily on the geometry of the system 


flow condition in such a “reactor” 
and on the rotational speed of the inner cylinder. 
Tayvon [2] has discovered that above a critical 
speed a regular pattern of toroidal vortices is set 
up which have about the same height as the 
distance between the inner and the outer cylinder. 
At much higher speeds the vortices become 
irregular, and a turbulence is 


obtained. Face [3] has shown that the axial flow 


more random 
velocity has only a secondary influence on these 
phenomena. In view of this available information, 
it could be hoped that an interpretation of 
measured diffusivities would be possible. 


Furthermore, the system chosen is of some 
interest for practical application mainly because 
a fairly uniform degree of agitation can be 
adjusted independently of the main flow rate. 
This special feature has been used, for instance, 
for liquid-liquid contacting such as extraction [4] 
and heterogeneous liquid reactions. An example 
of a reaction in this kind of equipment has been 


given by Revauurer [5] for the continuous 
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nitration of mono-chloro-mono-nitrobenzene with X ; 

7 s xe to (4) 
mixed acid. Paniku [6] has reported measure- X, 

nts on the heat transfer between the flow of  ,, ’ 

paceman: . hus, in order to compare «a and ¢ with the theory, 
water and oil and the outer cvlindrical wall in a , > aie : 
: ; og the complex transfer function X,/X, has to be 
similar apparatus. For such applications it ts 
desirable to know whether or not the longitudinal 
mixing impairs the effect of the continuous 


calculated from (3). 
The appropriate boundary conditions are : 
operation. XxX Xy 
CALCULATIONS (5) 
Suppose we have a reactor through which a liquid 
flows with uniform velo: wt) ‘a (Fig 1). The flow The latter condition has been discussed by 
is such that a constant diffusion coetlicient, D, PES “ ' = : : 
Daxckwerts [7] and also in [1]. The solution of 
exists in longitudinal direction. If we have the (3) and (5) is 
liquid flowing in with a variable concentration, C, , 
2 


the following differential equation obtains : > p er 9) oa 
“ 1 


\ . 
o¢ “ ; P2 e'Ps- pL 
oY Pi 


-~uU 


MASA KAP 
RRR RK 


pie 
KKK 








Fic. 1. One-dimensional flow with longitudinal diffusion. 


Here, C can also be considered to be the deviation with 


from the average value of the concentration. If - —- 
Pi & 4 y iwD _ 
2D ; 


C varies sinusoidally, we may substitute : ' 
P2 uz 


C= X=". 2 : , 
(2) It is easily shown that : 

thus representing C by a vector X rotating in the ‘ 

complex plane with an angular frequency w. With Lim —* "a, ora =1 and ¢ 


(2) equation (1) becomes an ordinary linear differ- 2° ““® 


ential equation : , : 
which represents a pure distance-velocity lag. 


po - ao 10x 0. ) A”: . 
dy 4 Lim —* = 1, ora = 1 and ¢ = 0, 


In order to determine D, the amplitude ratio, p+e Xo 


a, of C between y = a and y = 0 is measured, which means that in a perfect mixer the con- 
together with the phase shift, ¢, between the two centration has the same value at every point. 

signals, at various frequencies. The values of the The general result (6) is somewhat complicated, 
vector X at y = a, andy = 0 are related to each but under certain conditions it can be simplified 


other via « and ¢ in the following manner : to : 
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(8) 


In the present work we chose a ~ $1, whereas 
ul, 2D was of the order of unity. In that case it 
can be shown that (8) is a sufficiently accurate 
approximation of (7) if the measuring frequency, 
w, is chosen to be greater than 3u? D. 

It is to be noted that eq. (8) is the solution of 
(8) if the second boundary condition in (5) is 
replaced by 

y= oo, X =O, 


Thus the finite length of the reactor introduces an 
end effect which can be circumvented by taking 
a sufficiently small with respect to L and by 
measuring at relatively high frequencies. 


The amplitude ratio, «, and the phase shift can 
be calculated from eq. (8) to be : 


a 


, 1 tur? 1/2 | 
x N . 1 aoe =f —] 
exp[ Nol vala/ Nett) 4} 


— bw? 7? fe 
¢ No J5(/1+ Np? -a 


where 


(9) 


ua a 
and +r=>-- 


10 
2D u (10) 


Np 
Fig. 2 gives as an example the dependency of 
x and ¢ on N,, for several values of wr. An 
experiment under a given set of conditions gives, 
for a chosen value of wr, a result for « and ¢. 
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Fic. 2. Calculated relations between a and Np and between ¢ and Np for different values of wr. 
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From either of these a value of Np can be found, 
(Np), and (Np), respectively. If (Np), and 
(Np), agree with each other, and if the same 
values of these quantities are found at other 
measuring frequencies, it may be concluded that 
the longitudinal mixing can be described as a 
diffusion process. From the lines in Fig. 2 it can 
be seen that generally the phase shift will yield 
more dependable values of Np than the amplitude 
ratio. This is especially true for relatively low 
frequencies where the curve for « shows a weak 
minimum. 


Water was fed to the reactor at constant flow 
rate and measured with a rotameter. With this 
main flow a relatively small pulsating stream of 
a concentrated NH,Cl-solution was mixed, so 
that a practically constant feed rate to the reactor 
was obtained with a sinusoidally varying salt 
concentration. The pulsating flow of the con- 
centrated salt solution was obtained by displacing 
the with oil, the flow of which was 
modulated by a combination of a small constant 
rate feed pump and a buffer volume closed at one 
end with a reciprocating piston. 


solution 


The frequency 





series: 1 
Symbol: o 





- 


| 








Fic. 3. 


Pant 

The experimental arrangement is shown in Fig. 3. 
The (height 38 cm) 
stationary outer cylinder of * perspex.” Two of 
these were available with inner radii of R, — 4-52 
and 6-35em. The rotating solid evlinder had a 
steel core ; the outer diameter, 2R,, was adjusted 
to different values by covering the core with a 
layer of modelling wax after which it was ma- 
chined on a lathe to the wanted diameter. The 
rotational speed of the rotor was adjusted by 
means of a variable speed drive. 


EXPERIMENTAI 


“reactor” consisted of a 


n' (Ra-RFR, 
2* p) v 


Experimental arrangement. 


and the amplitude of the piston could be adjusted. 
The salt concentration at the inlet of the reactor 
was purely sinusoidal. The same applies to the 
electrical conductivity of the solution which at 
the concentrations used (up to 1-5 g/l) is linear 
with the concentration. 

Near the inlet and about half way the reactor 
a small conductivity cell was attached to the 
outer cylinder. The distance between the cells, a, 
was 16-3cm. The conductivity of the liquid in 
either cell could be recorded. By comparing the 
records of both cells under a given set of conditions 
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Table 1. Summary of results for water at room temperature. Measuring distance, a, between the cells 16-3 cm. 





| Dg D, D 
(10 on) (10 4 m? D n’ (Rg — R,)R, 
sec 


Series Ry Ry n (10 4 
- R,)R,| 4(Re + Ry) v 


m? 
(m) | (m) »/min) 


| 
sec —} n’ (Ry 


= —— ——___—__—_}—_—_—_ —_ a 


22 0-3 0-3 0-16 | 140 

70 0-85 | 0-9 0-15 440 

135 21 1-9 0-17 | 870 

0-0452 | 0-0249 250 2-3 2-4 0-11 1610 
390 2-6 3: 2-8 0-085 

520 3-1 | & | «ll 0-095 

735 42 5: 48 


0-3 ! 0-3 
10 1-2 
1-6 
21 

0-0452 | 0-0157 ‘ | 29 
41 
4-2 : 

48 5 5” . 2380 








Vl ! 120 
220 

3-1 ‘ 350 

| 00077 Hy 3-4 3 460 
40 li 570 

6-9 








| 00452 | 0-0051 





0-04382 





0-034) 


0-0234 


135 
0-0452 | 0-0157 250 
450 





48 x 10° m?, sec tv =055 x 10 ® m? /sec 


the amplitude ratio zand the phase shift ¢ between conductivity measurement in itself had a precision 
the two signals were measured. Although the of about 1%, the error in a, as derived from the 
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recordings, is estimated to be 5°%, or even greater 
for high damping ratios. The accuracy of ¢ was 
somewhat better but it decreased with higher 
measuring frequencies. The error in wr was not 
greater than 1%. 

An experimental run under a given condition 
consisted of the determination of « and ¢ for 8 or 
4 values of wr chosen between 1 and 10. By then 
comparing the values found for « and ¢ with the 
calculated graphs of « and ¢ vs. wr. for different 
values of Np an estimate was obtained for the 
values of (Np), and (Np), respectively. As can 
be seen from Fig. 2 a variation in ¢ of, say, 1% 
may cause a variation in (Np), between 3 and 
The 


determination of (Np), from a values is still more 


10%, depending on the value of wr used. 
uncertain. Consequently, the possible error in 


D, may sometimes be about 15% and in D, even 
as high as 30°. 

In Table 1 all results of these experiments have 
been assembled. The D, and D, values have been 
given separately. Except for a few cases the 
agreement between D, and D, is fair, considering 
the uncertainties This in general 
justifies the description of the longitudinal mixing 


by means of an effective diffusivity. For further 


involved. 


calculations an average value D, between D, and 
D, has been used with preference towards the 
value of D,. 

The linear velocity of the upward flow of the 
liquid varied in these runs between 0-002 and 
0-007 m/sec, which is small compared with the 
circumferential rotor speeds used (0-05-6 m/sec), 
and it could be expected that in these cases D 
would not depend on u. In 10 separate experi- 
ments under widely different conditions where 
besides the frequency also the main flow velocity 
was varied, no significant influence of & on D 
could be detected. 


Discussion 
All measurements were taken at a higher rotational 
speed than the critical value at which the vortex 
pattern appears. For this critical speed, n, 
60 n,’, TayLor’s formula [2] gives as a_ first 
approximation ; 


¢ * — P 6-5 
" . 5. 
v ( R, + R.) 


(11) 


For the present experiments the calculated n, 
varied between 0-5 and 6 rev/min, whereas the 
lowest value of n used was 22 rev/min. 

From visual observation with water to which 
some aluminium powder had been added it was 
ascertained that the height of the vortices was 
somewhat greater than the width of the annulus 
(R, — R,) only slightly with 
increasing rotational speed. Such a change in 
height of the vortices was discontinuous, since 


and increased 


the reactor can only contain an even number of 
At the same time it was found that the 
direction of the striae to be observed near the 
outer wall did not change much with the speed 
of the rotor. This indicates that, in first approxi- 


them. 


mation, the circumferential speeds of the rotor 
and of the vortices are proportional to each other. 
In order to relate the apparent diffusivity to 
the flow conditions this transport coefficient has 
to be compared to the product of a characteristic 
length and a significant velocity, for which then 
should be taken the height of the vortices and 
Instead of 
these quantities the product of (R, — R,) and 
of n'R, was chosen to form the dimensionless 
D/n'R, (R, — R,). As can be seen from 
Table 1 and from Fig. 4 this ratio is approximately 


their speed of rotation respectively. 


ratio : 


constant over the lower range of rotor speeds 
used. In this range 70% of the values of 
Dn'R, (R, — R,) lie between 0-14 and 0-19. The 
large deviations must be attributed to random 
errors, since no systematic trend could be de- 
tected. 

Two runs in which the viscosity was increased 
5-fold by using a glycerol-water mixture (Table 1, 
series 8) show no appreciable deviation from the 
other results. This indicates that the transport 
is entirely convective and that the molecular 
diffusivity of the salt molecules in the liquid has 
no influence once the toroidal flow pattern exists. 
The reported values of D are higher than the 
molecular diffusivity by a factor of about 10°. 

When the rotor speed is increased a transition 
can be observed from the regular vortex pattern 
to a more general turbulence. For a few cases 
the rotational speed, n, at which this gradual 
transition occurred was estimated visually. It 
can be seen from Table 2 that the dimensionless 
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group in TayLor’s criterion, eq. (11), has no 
significance here and that the width of the 
annulus apparently has a relatively stronger 


influence on the onset of general turbulence. 


Table 2. Rotor speed, n,, at which general turbulence 


sets in; water, room temperature. 





00452 00452 00-0452 


| 
—— eee 


0-0249 | 0-0157 0-0077 


Rg, m 


Ry _ R,, m 


150 


300 > 900 


. ! 
n,, rev/min | 
| 


. f 2 
n, (Ry— R,)R, 2%—™)' ‘188 


nm, (Ry — Ri)R, /, Ra- =} 
¥ Ry + R, 





Table 2 


general turbulence to appear between a rotating 


suggests the following criterion for 


inner cylinder and a stationary outer cylinder: 


period of 
-oscillation 
+$-925 sec 
—__ >< 


- 








salt sol'n 


2n,(R, eR, ~ 1.000. 


(12) 
(Re, r R,) 


In Fig. 4 the dimensionless diffusion number 
has been plotted against 2n’(R, — R,)?R, v 
(R,—R,). It is that the 
D/n'(R, — R,)R, significantly decreases as the 
criterion of eq. (12) is passed. This is in accord- 
ance the and the 
relative lowering of D may be mainly attributed 
to the decrease of the scale of turbulence in this 


evident value of 


with visual observations 


region as the rotor speed is increased. 

The above conclusions are far from being exact, 
but they may provide a basis for estimating the 
importance of the back-mixing effect if a physical 
operation or a chemical reaction is carried out in 
this kind of apparatus. Apart from the character- 
istics of the process involved, the value of 
uL,/2D (L = length of reactor) indicates whether 
the influence of longitudinal diffusion need to be 
If uL/2D > 1, e.g. 10 or more, the 
influence of back-mixing on the result of the 


considered. 


20-900 rpm 


> aa cells 


o--frer--7} 






































water 


pq 28-50 i/h Pm ’ 


to drain recorder 


Fic. 4. Proposed general correlation between the apparent diffusivity D, the 
dimensions of the annulus and the rotational speed of the inner cylinder. 
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process will generally be negligible. This applies 
to all of Partkn’s heat transfer measurements [6]. 
If on the other hand, for instance, with slow 
reactions and long residence times, 
then the 
In such a case back-mixing can 
mounting 


chemical 
reactor behaves like a 


uwL,2D <1 
perfect mixer. 
be suppressed to some extent by 
horizontal rings either on the rotor or on the 
stator which just leave sufficient clearance for 


the main liquid stream to pass but which prevent 
back-mixing between the different compartments, 
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NOTATION 


a = distance between points of measurement (m) 
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Flow of granular material through a circular orifice 
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Summary—Data have been obtained for the flow of granular materials through horizontal 
circular orifices, and correlated in the form of an empirical equation and nomograph with a 


70 


mean deviation of + 7 Data for inclined circular orifices are also represented by the same 
expression, when modified by a term involving the angle of inclination of the orifice as the only 
additional variable. The variables of importance influencing mass flow rate are the orifice diameter, 
particle diameter, particle density, and the angle of repose of the material, provided bed height and 
column-to-orifice diameter ratio exceed certain minimum values. 

Particle size ranged from 0-03 to 0-2 inches, particle density from 7-3 to 676 Ib./ft.*, and 
orifice diameter from 0-236 to 2-28 inches. 

Previous results for the flow of granular materials through restrictions are reviewed, and are 
found to be in reasonable agreement with the present correlation. 


Résumé 
et met les résultats sous la forme d'une équation empirique et d'un abaque avec une précision 
de + 7%. 


Les auteurs étudient lecoulement de matériaux granulaires au travers d' orifices circulaires 


Pour des orifices circulaires inclinés, la méme expression reste valable en tenant 
compte de cette inclinaison par une seule variable additionnelle. Pourvu que le lit ait une 
certaine épaisseur et que le rapport du diamétre de la colonne a celui de l'orifice dépasse certaines 


valeurs minimum, les variables qui ont une influence importante sur la vitesse d’écoulement sont 


le diamétre de lorifice, le diamétre et la densité des particules, langle de repos du matériau. 


Ces variables sont limitées aux valeurs suivantes : 


dimension des particules 0,03 a@ 0,02 inch 
7,3 a 6,76 Ib/te® 
0,236 a 2,28 inch 


densité des particules 
diamétre de lorifice 


Les auteurs trouvent des résultats en accord 


The need for quantitative information on the 
variables determining the rate of flow through 
a restriction has grown with the advent of large- 
scale processes employing moving beds of granular 
catalyst or other solid. Knowledge of fluid flow 
principles is of little help in such problems, 
since granular material systems, cither static or 
moving, deviate markedly from such principles. 
This is due mainly to the bridging of solid resulting 
from fhe forces exerted by a particle on neigh- 
bouring wall or particles. 

Previous investigations of granular solids flow 
have been reported. The least generalized of 
these was that of Newron, Dunuam, and Simpson 
[1] on the flowrate of 0-1 to 0-2 inch cylindrical 
catalyst particles 


clay cracking through a 


horizontal orifice of a flat bottomed container. 


(0,076 a 00,0508 em) 
(0,117 a 0,107 g,em*) 


(0,6 a 35,8 em) 


avec les données antérieures. 


The angle of repose of the material was reported 
as 88° and the minimum orifice diameter about 
six times the particle diameter. The empirical 
equation proposed relates mass flow rate, VW, 
to orifice diameter, D,, and bed height, H, as 
follows : 
iV = 8-50 D2 H°* (1) 
Since this equation is specific, it gives no informa- 
tion about the effect of particle diameter, angle 
of repose, or other possible fundamental variables. 
DeminGc and Menrine [2] studied the flowrate 
of a variety of materials through an inverted 
truncated cone orifice. Flow rate was found to 
vary with a power of the orifice size, and to be 
influenced by the size and apparent density of 
the particles, the angle of repose of the materials, 
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and the cone angle. The equation obtained would 
not be expected to apply to an orifice through a 
plate, that is, a cone angle of 180°. The equation 


was of the form : 


My 
D,?*D, 


[o-201 + (89-2 + 2-58 sin } ¢) 


(2s + 0-130 — 0-161 ms) (2) 
D, 
which relates time of flow, ¢, to the tangent of 
the static angle of repose, p,, orifice diameter, 
D,, particle diameter, D,, and the cone angle, ¢. 

The most generalized investigation was re- 
ported by Takanasi [3] who studied many types 
of sand, shot, vegetable seeds, and mixtures of 


His 


result was expressed by the empirical equation : 


- D, : 
= fem [seore(n)] 


Where a and b are constants dependent only 
on the system of units used: g is gravitational 


sand and sulphur and sand and talcum. 


acceleration, f(y.) is a power function of the 
tangent of the surface kinetic angle of repose, 
and other terms are as previously defined. 

The appearance of the D, in two terms tends 
to produce a curved line on a log-log plot of flow 
rate vs. D,. Some examples of the observed data 
given appear to have some curvature in such 
a plot, but most could equally well be expressed 
as a Straight line, especially when one considers 
the size of the apparatus used and the resulting 
probable accuracy. It is also quite likely that the 
data on large particles were excessively influenced 
by the container walls, which Takanasi did not 
consider. In addition there is some uncertainty 
whether reciprocal volume flow rate is to be 
based upon bulk volume or particle volume. 

In view of this meagre and uncertain informa- 
tion, the present investigation was begun to study 
the effect of the following variables : (1) orifice 
diameter, (2) particle diameter, (3) particle 
density, (4) angle of repose, both static and 
kinetic, (5) bed height, (6) void fraction, and 
(7) angle of inclination of the orifice. Preliminary 
study had indicated flow rate could also be 


FRANKLIN and L,. 


N. JOHANSON 


influenced by the following factors : (1) particle 
to orifice ratio exceeding 0-2, (2) orifice to column 
diameter ratio exceeding 0-16, and (3) bed height 
less than one column diameter. Correlation was 
therefore restricted to values within these limits. 
In addition, particles selected were of sufficient 
size that effects of electrostatic 
attraction or repulsion were negligible. All par- 


ticles were air dried to avoid agglomeration. 


and weight 


EXPERIMENTAL 


Matcrial was poured loosely into vertical glass 
columns and allowed to flow downward through 
a circular orifice in a ¥ inch plate. The orifice 
was horizontal and centered at the column bottom 
except for the runs investigating the effect of 
angle of inclination from the horizontal. Flow 
rate was obtained by determining the weight 
of material flowing in a given time interval, 
measured by spotwatch. This method allowed 
several rate values to be obtained during a run 
as a test of rate variation with decreasing bed 
height. Column diameter ranged from 1-76 to 
8-81 inches depending upon the diameter of the 
orifice used. Four smaller columns were used to 
study the wall effect at constant orifice diameter. 
The inner surfaces of the columns were roughened 
with drafting tape to study the effect of wall 
roughness. 

Particle diameter ranged from 0-03 to 0-2 
inches. The diameter was determined by the 
average of direct measurement for large particles, 
and by arithmetic average screen mesh size of 
the fraction for smaller particles. 

Particle density was measured by the immersion 
method, the fluid used depending on the porosity 
of the particle. It was felt that actual particle 
density, including internal pores was a more 
fundamental variable than either bulk or true 
material density. Bulk density and the related 
void fraction were found to be dependent on the 
manner of placing the material in the column, 
whereas the flowrate, once started, is not so 
dependent. It has been found also that the void 
fraction of the static bed is not necessarily the 
same as for the moving bed. On the other hand, 
true material density could not be correlated 
with flowrate variation. This appears to indicate 
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that the flowrate is a volumetric, particle-spatia] 
phenomenon. 

For porous materials mercury was used as the 
immersion fluid, with a vacuum being utilized 
to remove air from the inter-particle voids in 
order to reduce the possibility of unfilled voids. 
For the non-porous particles water was used. 
A thin wax coat of negligible volume was applied 
the 


mercury 


to the puffed rice making possible 
water as fluid since would 
both the 


The particle density was varied from 


Immersion 


have penetrated and compressed 


particles. 
676 |b. per cu. ft. for lead shot to 7-34 Ib. per cu. ft. 
for puffed rice. 

Static bed void fractions were calculated from 
particle density and bulk density values for the 
static bed. 

The angle of repose would appear to be an 
important characterization factor for less easily 
defined particle properties, such as shape, rough 
effective and effective void fraction. 


ness, SIZC 


Previous workers have used the static angle of 


kinetic 


may be 


the surface 


first 


and 


The 


repose [2], 


repose [3]. defined as the 


maximum angular deviation from horizontal of 


the surface of a static pile of granular material 
The 


is the slope of a continuously moving surface of a 


Its tangent is defined here as jp,. second 


pile of granular material 


(tangent fg). 


Observation of the “shear angk hetween 


relatively fast moving end slow moving particles 


as granular material cones toward an _ orifice 


suggests that such an internal angle may be more 


proper fundamentally. Though this angle was 


not measureable with any precision in the columns 


used, an internal angle’ believed — similar 


(tangent = p,) was easily and reproducibly 


measured in an apparatus [4] constructed for this 
This used for 


determining static and surface kinetic angles of 


purpose. apparatus was also 
repose. 

The apparatus consisted of a rotating drum one 
foot in diameter and eight inches long, with axis 
fixed in a horizontal position, and having a 
front face of transparent lucite plastic through 
which the various angles of repose could be 


observed. The curved surface was roughened to 


prevent slippage of the material along this wall. 


use of 


angle of 


The drum could be rotated at constant speeds, 
from very low rates to about five or six r.p.m. 
The maximum angle of repose attainable by the 
solid before moving or cascading downward was 
Kinetic 


angles were measured with the particles rolling 


recorded as the static angle of repose. 


downward over the slanting surface at essentially 
a constant rate. The internal kinetic angle with 
tangent defined here as yp, could be observed as 
the 


moving stream and the stream being carried 


the plane of contact between downward 


upward by the motion of the rotating drum as 


shown in Fig. 1. The use of internal kinetic 





Surfoce Formed 
by Downword Moving 
Porticies —2 


Surface Kinetic internal Kinetic 


Horizontal 


Plone ot Which 
Particles Change 


Arrows Denote 
Pattern of Particle 
Motion 


Drum Rotation Direction 











Fic. 1. Surface kinetic and internal kinetic angles of 


repose. 


angle of repose gave a better correlation of the 
flow rate data than did use of either static or 
surface kinetic angle. This may have been due 
in part to the dependence of surface angle upon 
the rate of rotation of the drum. The effect of 
the lucite face on the angle of repose appeared 
negligible, since the surface at the lucite face 
was not appreciably deformed from that of the 
bulk of the material in the drum. Intercorrelations 
of the three angles of repose; static, surface 
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kinetic, and internal kinetic are shown in Fig. 2. 
Although this correlation is of adequate precision 
for this work, it should be accepted with caution 
in view of the limited knowledge of the factors 
influencing the various angles of repose. 
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Fic. 2. Inter-relationship of angles of repose 


Results obtained for cach of the materials used 
are summarized in Table I for horizontal! orifices, 
and Table II for inclined orifices. These will be 


discussed separately. 


HoriIzONTAL OrtFick RESULTS 


Preliminary experiments were conducted to 


determine the influence of the nearness of the 
column wall upon flow rates. Results obtained 
varying column size, for an orifice to particle size 
consistent flow, 


minimum for 


There is an influence on 


ratio near the 
are shown in Fig. 3. 
flow rate as orifice and column diameters approach 
the same value. It is believed that the difference 
in radii of orifice and column in terms of particle 
diameters is the pertinent variable. Thus the 
abcissa expressed as (D. — D,)/D, may be a 
more general expression of this effect. This term 
and the one Fig. 3 are 
(D, — D,)/ D, D,/D,. 


related ; 


When 


used in 


(D,/D, — 1) 


JOHANSON 


the difference between orifice and column dia- 
meter is greater than about 30 particle diameters, 
the influence of the wall is negligible. In most 
applications it is relatively easy to obtain this 
condition, and this was done for all subsequent 
measurements in this work. This relationship 
between column and orifice diameter was only 
slightly influenced by the roughness of the wall. 
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Fic. 3. Influence of container wall upon flowrate. 
Straight lines were obtained when the logarithm 
of the flow rate was plotted versus the logarithm 
of orifice diameter, so that the influence of the 
orifice diameter is completely satisfied by one 
This is shown in Fig. 4. The 
slopes of the lines and hence the exponents of 
the diameter varied the material 
used from 28 to 3-1 with an average 
value of 2-93. A plot of the exponent versus the 


simple factor. 
orifice with 
about 


static bed void fraction indicated that there was 
some degree of correlation, the exponent tending 
to increase with void fraction, but the points 
were scattered sulliciently that an attempt to 
express an exact relationship was not believed 
warranted, 

Having determined the influence of the orifice 
diameter, the other 
variables would appear in the intercept values 
at an orifice diameter value of unity in such a 
log-log plot. Expressing these intercept values 
of flow rate in units of particle volume rate no 
further correlation could be found between these 


remaining influence of 


values and particle density, void fraction, or 
bed height, provided that the bed height was 
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Table I. Horizontal orifice results. 





Material | ‘* | W (Bq. 5) 





Glass Beads, 16-28 mesh 0-236 
(Scotchlite) | ©2907 
(422 
0-614 
1-051 
Glass Beads, 5 mm borosilicate | 1°35 | 0-196 
2:28 
Ion Exchange Resin spheres, 0-236 0-031 92 21°5 0-307 
16-28 mesh 0422 | 1-69 
0-799 10°27 
Lead Shot, ** Dust” 0-236 00-0425 | 2) 1-85 
0-207 3-83 
0-422 11-92 
0-800 72-1 
Lead Shot, No. 12 0-422 0-050 114) 
| 0614 32-1 
0-800 
Cracking Catalyst, 
Used Bead Spheres 
Through 6 mesh 0-799 0-398 
0-925 
O51 
‘378 
4-6 mesh ‘378 
“693 
“930 
Unsized ‘378 0-130 0-398 
“930 
Cracking Catalyst, 
Clay Cylinders 
New Unsized ‘B52 
‘28 
Used Unsized ‘378 
‘93 
‘28 | 
Puffed Rice ‘35 0-205 
93 
-28 
Crushed Olivine Rock 
16-28 mesh 0-236 0-031 | 0-547 
0-422 
0-800 
4-6 Mesh | 1-35 0-158 0-535 
2-28 | 
Coal (Tono), 16-28 Mesh | 0-236 0-031 0-530 
0-422 
0-614 
0-800 





| 





* Each value of W represents the mean of two or more determinations, with an average deviation of less than two 
per cent. 
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Table II. Inclined orifice results and comparison with equation 6. 





Flow Rate, lb./min. 
Experimental Results 
@=30° 06=60° @=9%90 


Material D, (in.) | Dy{in.)| 8 =0 


4°91 
2-39 
0-717 


0-50 0-031 
0-375 
0-25 





Glass Beads, 
16-28 mesh 
“ Scotchlite * 
0-50 0-0425 
0-375 
0°25 
| 


Lead Shot, “ Dust ” | 


4°47 


0-634 


Flow Rate, lb./min. 
From Equation 6* 
@ = 30° @ = 60° g =90° 


2-33 
1-12 
0-340 


8-62 
1-76 
0-529 


2-01 
0-873 
0-254 


3-58 


0-440 0-667 

10°3 
4:12 
1-24 


15°85 
6-35 
1-90 


19-9 
7-96 
2-39 


11-23 
3-83 
1-21 


16-85 





* Based upon data of column 4 for @ = 0°. 


greater than the height of the cone formed by 
the moving particles as they approached the 
orifice. 

No influence of excessive compaction, bridging, 
or related effects were observed as bed depths 
were increased, the six foot 
depth studied. This is in agreement with previous 
work [1, 3], and may indicate that the wall 
supports the major weight of granular material, 
eat least in the relatively narrow columns usd in 


up to maximum 


such work. 

In a_ preliminary the 
reciprocal of the intercept values obtained from 
Fig. 4 versus particle diameter, it was found that 
the data could be represented by straight lines 


correlation plot of 


if the tangent of the internal kinetic angle of 


repose yz, Was used as a parameter, in accordance 
with the following equation ; 


(; B) /(D, A) =C yp, + E. (4) 


Depending on the values of A and B the family 
could either be a series of parallel lines or could 
actually intersect at a finite point. Upon re- 


arrangement the equation may be expressed ; 


1 
kK 


or y = az, + ba, + cx, + d. 


= Cys D, + ED, —C A pg —(A E—B). 


where y = 1/K, 2, = 3 D,, % = D,. @ = py, 
a=C,b= E,c = —CA,andd - AE + B. 
The coefficients a, b, c, and d were evaluated by 
the method of least squares and incorporated 


into the basic equation to yield the final empirical 
result : 


pp D,?* 


(6-288 py, + 23-16) (D, + 1-889) — 44-90 


If eq. (5) is plotted with D, as abcissa on log-log 
paper all the data can be represented by a single 
straight line with a slope of 2-93, the maximum 
12%. mean 


deviation being about and the 
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Fic. 4. Influence of orifice diameter upon flowrate. 
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Fic. 5. Proposed flowrate correlation. 


deviation 7%. This correlation is shown in Fig. 5 
and Table I. 
structed on the basis of the above equation to 


The nomograph, Fig. 6, was con- 


allow more rapid calculation. 


INCLINED ORIFICE RESULTS 
The effect of the angle of inclination of the orifice, 
as summarized in Table II, is from the work of 
GuNDERSON [5]. The lead shot and the glass 
beads used were the same as used for horizontal 
orifice studies, and the apparatus and method of 
determining flow rate were similar. The orifices 
used were approximately }, ¢ and 4 inch, drilled 
through »; inch galvanized sheet metal. The 
orifice plates were mounted at the truncated end 
of a two inch diameter cardboard tube, for 30° 
and 60° inclination. For 90° inclination the 
orifice was mounted in a side face of a rectangular 
box of the same diameter as the circular cylinder 
below which it was attached. 

As illustrated by Fig. 7 for lead shot, the data 
can be represented by straight lines if mass rate 
of flow is plotted versus the cosine of the angle 


W, 
(b./mia. 

















FROM OC AND Dp, K)*Kp IS OBTAINED 
PROCEED TO W KNOWING Do AND Mp o)_) 
OR TO Do KNOWING W AND Op 











Fic. 6. Nomograph of proposed flowrate correlation, 


equation (5). 


the orifice. 


An extrapolation of such lines to intercept the 


of inclination from horizontal of 
abcissa at zero flow rate suggests that flow rate 
will approach zero as the axis of the orifice 
approaches a line perpendicular to the kinetic 
angle of repose. Since this is intuitively logical 
the data were related to the flow rate through a 
horizontal orifice by means of a multiplying 
factor involving the angle of inclination of the 
orifice, and the angle of repose of the granular 
solid, thus 

w, = W, cos a + cos ] (6) 

cosa + 1 

where. W, is the flow rate through an orifice 
inclined at an angle @ with the horizontal, and 
a is the kinetic angle of repose. The factor 
reduces to one for @=0, and to zero for 
cos 8 = — cos a, thus satisfying the two extreme 
conditions. Fig. 8 and Table II compare the 
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Fic. 7. Flowrate correlation of lead shot through inclined 
orifices. 


1.0 


experimental data for inclined orifices with eq. (6). 
The equation appears to fit the data within the 
experimental error. 

Thus the flow rate through an inclined orifice 
can be estimated from a knowledge of the flow 
rate through a horizontal orifice, and the inclina- 
tion angle of the orifice without the need of 
additional empirical constants. Flow rate through 
the horizontal orifice can be experimental, or 
obtained by using Equation 5 if the 
diameter and particle properties are known with 
for the 


orifice 
sufficient accuracy. GUNDERSON’s data 
flow rate through a horizontal orifice were used 
as the basis of comparison in Fig. 8, since the use 
of the nominal values of orifice diameter together 
with eq. (5) resulted in values for W, which were 
10% to 20% 
results. 


from lower than experimental 

It would be of interest to compare the results 
of previous investigators with those predicted 
by means of eqs. (5) and (6). No data other than 
GUNDERSON’s is available for through 
comparison of 


flow 
inclined orifices and hence no 


eq. (6) with other workers’ results is possible. 


N. JOHANSON 
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Fic. 8. Proposed flowrate correlation for inclined orifices 


Equation (5) for horizontal orifices can be com- 
paried with previous work, however, and the 
results of such comparison are summarized in 
Table IIL. 


comparisons are discussed below 


The assumptions involved in these 


TAKAHASI does not values obtained 


experimentally, except in the form of plotted 
points for obtained 
from his Fig. 7 involving a thousand-fold variation 


report 


five materials. Values 
in flow rate are of questionable accuracy, but 
have been compared with those predicted by 
his equation and by eq. (5) for the materials 
quartz sand, “iron sand,’ and lead shot. In 
lieu of a reliable particle density for “ iron sand,” 
a void fraction of 0-46 was assumed, and particle 
density calculated from bulk density. Internal 
kinetic angles of repose were obtained from 
Takahasi’s surface angles by means of Fig. 2, 
for quartz and iron sand. The value for lead 
shot is from this work, Table I. 

The agreement between experimental results 
and either equation is less than would be desired, 
but is suflicient to indicate eq. (5) represents 
the about TAKAHASI'S 
equation. Much of the discrepancy can be 
attributed to inconsistencies in the experimental 
data. In using these data, care was taken to 
avoid those results for which the ratio of particle 


data as well as does 
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size to orifice size, or orifice size to column size, 
was too great to assume independence of flow rate 
from these factors. 

It is believed that the clay catalyst upon which 
eq. (1) is based is substantially the same as the 
new clay catalyst of this work. Properties of the 
latter catalyst were therefore used to enable a 
comparison of eq. (1) and eq. (5) with experimental 
results of the present work. Use of eq. (1) with 
bed height equal to one foot results in values 
about 17% high, while eq. (5) results in values 
about 8% high. In view of the assumptions 
necessary to use eq. (1), the two equation may be 
considered to be in agreement. 

DemiInc and MEHRING recommend their equa- 
tion only for cone orifices, of cone angle ranging 
from 20° to 110°. Their data were obtained for 
30°, 60° and 90 
flow rate with increasing cone angle is not marked, 
and it would be expected that application of their 


equation to cone angles approaching 180°, or a 


cone angles. The decrease in 


horizontal orifice, would not be greatly in error. 
The 
angles of 180 


was actually made for cone 
2a, where « is the kinetic angle 
the 


material 


comparison 
of repose. This represents approximately 


natural formed when granular 
issues through a horizontal orifice. It is reasonable 


to expect that flow rate would not be influenced 


cone 


by cone angles greater than this, and therefore 
eq. (2) which is that used by Deminc and Menrine 
need be extrapolated only to this point for proper 
comparison with eq. (5). 

There is some difficulty in making a direct 
comparison between their results and the present 


equation even if such extrapolation ts justified, 


however. Their equation is based upon static 
angle of repose and bulk density of the granular 
solid, and they did not report values of particle 
density or fraction void and the kinetic angle of 
Values from Table I 
were therefore used for glass beads and lead shot. 
Kale and mustard seed densities were obtained 
from bulk density and an estimated void fraction 
of 0-4, and published densities for urea and 
Internal kinetic 


repose required in eq. (5). 


phosphate rock were used. 
angle of repose was estimated for these materials 
With the exception of 


which is considerably 


by means of Fig. 2. 
Phosphate rock No. 2, 


FRANKLIN and L. N. 


JOHANSON 


finer than any material used in the present work, 
Table III indicates the equations are in reasonable 
agreement at the point where they are each 
applicable, namely for cone angles of 180 2a. 
The values predicted from Deminc and Meurina’s 
equation are also in line with their experimental 


results for 30, 60 and 90 degree cone angles. 


CONCLUSIONS 


The mass rate of flow of granular material through 
circular orifices in thin flat plates can be calculated 
from a knowledge of a few easily measurable 
properties of the system. These are the orifice 
and particle diameters, the particle density, the 
angle of inclination of the orifice with the hori- 
zontal, and a kinetic angle of repose of the 
granular material. A device is suggested for 
measuring the latter value. Equation (5), or the 


analogous nomograph, Fig. 6, can be used to 


calculate flow rate through a horizontal orifice, 
and eq. (6) can be used to modify this for inclined 
orifices. For cone-bottomed containers terminat- 
ing in a horizontal orifice, eq. (2) from the work 
and MEHRING used up to 


of Deminu can be 


cone angles of 180 2a, where o is the angle 
of repose. For cone angles greater than this, 


eq. (5) is satisfactory. 


NOTATION 


Diameter, inches. 
o to orifice and p to particle. 
Height of bed above orifice, feet 


D, 1 of flow rate 


Subscript c refers to column, 


Intercept at versus orifice 
diameter plot. 
Time of flow, minutes /unit volume. 


flow, ft.* 
through column 


min Subscript « 


rate of 
flow rate 


Volume 
refers to 
D,/ D, 
Mass rate of tlow, lb. 
script o refers to flow through horizontal orifice ; 
subscript @ to flow through inclined orifice. 
Angle of repose. Refers to internal kinetic angle 
unless otherwise stated. 

Void fraction of static packed bed. 

Angle of inclination of orifice from horizontal. 
Tangent of angle of repose. Subscript 1 refers 
to surface static angle ; 2 to surface kinetic angle, 
and 3 to internal kinetic angle. 

Density of material Ib./ft.°. Subscript 6 refers 
to bulk density ; p to particle density. 

Cone angle for conical orifice. 


for which 
> &®. 


min. No subscript or sub- 
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Dimensional systems and systems of units in physics with special 
reference to chemical engineering 


Part I. 


The principles according to which dimensional systems and systems 


of units are constructed 


A. KLINKENBERG 
N.V. De Bataafsche Petroleum Maatschappij (Royal Dutch/Shell Group), The Hague 


( Received 23 January 1955) 


Summary 


In Part I the principles are examined of building systems of dimensions and units 


in physics with special reference to the aspects for the chemical engineer. 
This part provides the background for Part II which gives : 


a formulation of practical rules, 


a discussion of the confusion in the field, and 
conclusions re advantages and disadvantages of various systems. 


Résumé 


Dans cette premitre partie, auteur étudie les principes rationels de l'établissmenet des 


systémes de dimensions et d'unités des grandeurs physiques au point de vue spécial du Genie 


CHIMIQUE. 


Cette étude servira de base & une seconde partie d'écrivant : 


(a) Vénoncé des régles pratiques. 


(b) la discussion des confusions existant présentement dans ce domaine. 
(c) Avantages et défauts des différents systémes actuels. 


1. INTRODUCTION 


In chemical engineering much of the knowledge 
which authors want to convey to readers is in the 
form of equations. It is believed that always the 
ultimate aim of such equations is their use in 
numerical calculations be it by the author or by 
the reader. 

It is the purpose of this paper to examine how 
this aim may be reached with the least risk of 
making mistakes. 

A well-known text book [1], page 29, says : 

“ The data of the chemical engineer will often be 
found in various systems of units and he must be 
able to disentangle them. The burden of straight 
thinking is upon him and he must be able to avoid 
the confusion induced by the terminology employed 
by brother engineers.” 


This quotation only covers the reader's aspect. 
The crux, however, lies with his “ brother chemi- 
cal engineer,”’ the author. While the philosophy 
of this paper is to allow the author the use of any 
system of units — disallowing of course any hap- 


hazard assembly of units - the point is equally 
strongly pressed that the reader has the right to 
see the system of units used correctly. 

This latter point carries added weight in view 
of the circumstance that the user of an equation 
often has no time to go through its complete 
derivation. In fact, he glances through the 
article, until he finds a combination of the right 
symbols, and he wants to use that equation 
rightaway, possibly with units of his own choice. 

For the above to be possible without undue 
effort it is necessary that dimensions and units 
should be well understood, handled correctly and 
indicated clearly. 

The following text gives in Part I a description 
of the principles used in building up systems of 
dimensions in current use. It has been tried to 
give this description with the least amount of 
arithmetic ; for this reason, some statements may 
appear to be without proof. A rigid background 
may be obtained with the aid of works by 
BripeMan [2], Focxen [3] and others. 
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In Part Il a number of practical rules are 
formulated for the use of authors and readers of 
papers and reports in which calculations of a 
These 
literature 


physical nature are carried out. are 


followed by a critical discussion of 
quotations. 

It is clear that the above aims can be reached 
much more efficiently if (1) symbols are standard- 
ized; (2) everyone uses the same dimensional 
formulae and units. The first of these points will 
not be touched. The second would have following 
obvious advantages: no attention need be 
devoted to mastering several systems; no con- 
version factors are required and no actual work 
of conversion need be done; a better feeling can 
be developed for the numerical values of physical 
magnitudes* if they are always expressed in the 
same units. 

In order to assist in the choice of the most 
suitable system a summary (Part IT, § 5) has been 
prepared of the advantages and disadvantages of 
various systems as derived from a prior’ consi- 
derations and as noted in the critical discussion 


of the literature. 


2. Quantities OCCURRING IN 
EQUATIONS 
The equations, seen as relations between numbers 
One way of looking at equations is that they 
represent certain mathematical relations between 
numbers. 
These 


namely : 


numbers may be of various types, 

(a) Simple numerical factors (definitional con- 
stants according to GUGGENHEIM [4)). 

(6) Numbers expressing the size of a physical 
magnitude and represented by a letter sym- 


bol. 


of the size of a physical magnitude to the size 


Such letter symbols stand for the ratio 


of the corresponding unit. 

Numbers expressing the size of a physical 
magnitude and given as such, i.e., as a 
numeric, the mentioned 
under (0). 


equal to ratio 


= Following Focken [3], length is a physical quantity, 
and the length of a particular object is a physical magni- 
tude. 


(d) So-called dimensional constants, represented 
by a letter symbol or by a number. 

These arise from the fact that the units for 
certain physical quantities can be derived in 
more ‘than one way from the units for the 
quantities assumed fundamental. 


A few simple cases will provide examples for 
the above four types of factors. 


1. Free fall without initial velocity 


If path is measured in metres and time in seconds, 
we may write : 


49 * (1) 


In this equation we find : 
(a) a numerical factor }; 


(b) sizes of physical magnitudes “ path” and 
“ time,”’ expressed by letter symbols s and ¢; 


(e) size of physical magnitude “ acceleration,” 
expressed by numerical value 9-8 ; 


no dimensional constant. 


2. Shear stress in a liquid 
The relation between shear stress 7 and velocity 
gradient du/dy in a fluid is represented in a 


certain case by : 


du 


&e 7 "7S 


when 7 is measured in |b.,/sq.ft.* 
» (= viscosity) in lb. /ft.sec.* 
u (= velocity) in ft./sec. 
y (= a coordinate) in ft. 


Besides sizes of physical magnitudes (see para. 
(b) above) the equation contains a dimensional 
constant g. of numerical value 32-2. 

The occurrence of this constant is due to the 
existence of alternate routes in deriving units. We 
might have defined the poundal (force acting on 
1 Ib. mass with an acceleration of 1 ft. sec.*) and 
used this unit in the equation without a factor 
g.. However, we the Ib. 
(= force acting on a 1 lb. mass with an accelera- 


have chosen force 


+ Ib.¢ =» pound force ; Ib.,, = pound mass. 
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Since the 


tion of g = 32-2 ft./sec.*) instead. 


number of Ibs. force is the number of 


82-2 
poundals, the factor g, = 32-2 
included. 


must then be 


3. The first law of thermodynamics 

If energy is expressed in two ways, viz. as work 
in mechanical units and as heat in heat units, the 
mechanical heat equivalent enters as a dimen- 
sional constant. 


Of the four types of numbers that may occur 
in equations, one, viz. type (c), is often undesirable. 
We may discuss this point by taking the equation 
(1), s = 49 as an example. This 
equation has the advantage of being immediately 


form of 


usable for anyone employing m and sec. as his 
fundamental units, } gt? the 
substitution of g On the 
other hand equation (1) has the great disadvantage 


whereas with s 
9-8 is still needed. 


of requiring conversion if it is to be used with s 
and ¢ expressed in units other than m and sec. 

For the generality of the treatment it is there- 
fore considered advisable to work as long as 
possible in letter symbols, to replace the constants 
such as g and various physical properties by their 
numerical values only immediately prior to the 
use of the equation for a numerical calculation 
and from that point onward, always to state 
precisely what units are to be used. 

Thus : 


s = 492 (units : m and sec.) 


exactly as one would say 
t = 20°C 
and never simply 
t = 20 


The equations containing nothing but purely 
mathematical numbers (a) and letter symbols 
representing the sizes of physical magnitudes (6) 
or the values of dimensional constants (d) are 
called complete. They are invariant for changes in 
the size of the units. Indeed, one should wish the 
description of the laws of nature to be independent 
of our choice of units. It is, however, to be noted 
that such invariance can be reached with regard 


to changes in size of units and fundamental units, 
but not with regard to changes in type (see 
examples 2 and 8). 


The equations, seen as relations between physical 
magnitudes 


So far, all symbols stood for numbers. One may 
also look at the symbols differently and consider 
them to represent the physical magnitudes them- 
selves. This, however, is the subject of violent 
attack by some authors, who maintain that it is 
impossible to multiply or divide physical magni- 
tudes. 

The stand taken here is that it is perfectly 
legitimate and very convenient to do so. 

This simply means that the concept of multi- 
plication is given a wider scope than it has in 
arithmetic. 

Even in pure mathematics this procedure is 


When making the 
b, each vector 0 is considered to be the 


known. vector product 
c=a 
product of a numerical value v and a unit vector 
v,. The equation 


cc, = aa, x bb, (8) 


then falls apart in the two equations 


(8a) 
and 


(3b) 


It might astonishing to derive two 
equations from one. This point will be clarified 
by an example from physics: 

Since a body in uniform motion at 10 m/sec. in 


1 hr. covers a path of 36 km one may say : 


1 hr. 


seem 


36 km = 10 m/sec. 
This statement is evidently highly impractical. 
It may be made practical in two ways, viz. by 
writing : 
Firstly : 
86000 m = 10 m/sec. x 8600 sec. 


Here the equation falls apart in two : 


36000 = 10 x 8600 
and 
1m = 1 m/sec. x 1 sec. 
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We have achieved this by relating the unit of 
velocity to the units of length and time. 


Secondly : 
36 km = 10 m/sec. x Lhr. x 3-6 umn.see. 
m.hr. 

The added factor is a dimensional constant. 

Again, both the numbers and the units are in 
order. 

In either case, by considering every physical 
magnitude to be the product of a numerical value 
and a unit, we have obtained the possibility to 
use the equation in two ways, namely to calculate 
numerical values and to carry out a check on 
dimensional homogeneity. 

The conversion into other units by varying the 
size of the fundamental units also becomes very 
simple, viz. : 

- 36 —____ — 
8600 sec. 

There is no risk of introducing conversion 
factors the wrong way round. 

The interpretation of the equations as relations 


between physical magnitudes is used or advocated, 
amongst others, by CorneLivs and Hamaker [5], 
FLEISCHMANN [6], GuGGENHEIM [4], Pace [7], 
Pout [8], and Tuorne and Watsnaw [9]. 

A difficulty arises if the equation would contain 


logarithms of physical magnitudes. On closer 
thought, it is seen only to contain logarithms of 
ratios, i.e. differences of logarithms. However, 
one such logarithm may be visible and the other 
concealed in an integration constant (examples : 
integrated Clapeyron equation; Arrhenius 
equation ; entropy of an ideal gas). 

This forces one, if one writes the equation in 
the usual form, either to accept logarithms of 
physical magnitudes and units or in this case to 
use numerical values only. 

The author is inclined not to take this difficulty 
too seriously and to argue that, whatever illogical 
or physically impossible is done, it is done twice 
and cancels out. Reference is made to BripGMAN 
[2] and Focken [3], page 43. 

Of the afore-mentioned two methods of obtain- 
ing dimensional homogeneity, the first by far is 
the most common, but it is absolutely essential 


for the understanding of systems of units that the 
second method is well recognized. 

So far, we made a direct jump to units, to be 
distinguished as fundamental units (in the first 
example : m and sec.) and derived units (in the 
first example : m_/sec.). 

In order to write and read complete equations 
in letter form, however, there need be no more 
than agreement on fundamental quantities, i.e. 
on types of fundamental units to be used (length 
and time in the first example), and agreement on 
certain conventions that fix numerical constants. 

Such conventions are, for example: deriving 
the unit of area from a square with unit side and 
not from a circle with unit diameter; deriving 
the Reynolds group from a pipe diameter and not 
from a radius ; choosing between placing a factor 
1/4m in the formula for the electrical capacity 
of a plane condenser (in the cgs system) or a 
factor 47 in that for the capacity of a sphere 
(“rationalization ’’ as applied in the practical 
system). 

The matter of sizes of fundamental units only 
arises at a later stage when numerical values are 
substituted and the resultant equations are no 
longer complete. 

These two steps should therefore in the deriva- 
tion of equations be made separately. 


CHOOSING 
UNITs 


PRINCIPLES OF 
OF FUNDAMENTAL 


3. SOME 
TYPEs 
The problem of units in physics and engineering 
is sometimes considered to involve primarily a 
choice between “the” metric and “the” English 
system. 
The present discussion will underline that this 
is a serious understatement. It should be noted : 


1. That the matter of choice of fundamental 

quantities (types of fundamental! units) is 
much deeper rooted and requires more in- 
sight than the matter of the choice of size of 
fundamental units (see, e.g., FLEISCHMANN 
(24}) ; 
That even for given fundamental quantities 
the metric as well as the English system allow 
various alternatives (em-m; g-kg; in.-ft. ; 
sec.-hr.) all of which are in use. 
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The above example of movement at uniform 
velocity showed two ways of obtaining a complete 
equation : 
Firstly : 

units of length and time. 


by relating the unit of velocity to the 


Secondly : by introducing a dimensional constant 
without relating the unit of velocity to the 
units of length and time. 


For this example every one agrees that the first 
method is to be preferred, since it is the simpler. 

In continuing in this manner to relate derived 
quantities to fundamental quantities one builds 
up a system of dimensions. 


The so-called dimensional formulae. such as 


LT’ 
[LT] 


velocity 


are nothing but a useful short-hand notation to 
see how this is done. 
The statement 


unit of velocity = m.sec. 


serves the same purpose. In view of § 2 this 
statement, however, at this point introduces an 
unnecessary limitation. 

It is preferred to call [L T''] the “ 
* dimension.” 


dimensional 
formulae ” and not the The latter 
term should be reserved for the set of exponents. 
Velocity thus has the dimensions +1, 0 and 
— 1 in length, mass and time. In this way a 


dimensionless quantity has dimensions zero 
in all fundamental! quantities, but its dimensional 
formula is unity. 

This principle of avoiding dimensional constants 
when carried through further and further, all 
through physics, finally defeats itself. 

It may be argued that one has gone too far 
already in the egs system. 

In this electric 
quantities are derived from length, mass and time, 


with the consequent appearance of weird dimen- 


system all and magnetic 


sional formulae such as 


electric charge in electrostatic system 
(I 3/2 Mi‘? te 
Such dimensional formulae do not, as it is 
sometimes believed, betray the deepest secrets 


of nature. Instead they are rather confusing 


and they have been an obstacle to generations of 
physicists. It will be seen that a system with 4 
fundamental units is more convenient. 

The argument further to 
extremes and in order to bring out very clearly 
the arbitrariness of all systems of dimensions, 


this will, in the following, be done. 


may be carried 


The starting point is that the cgs electric and 
magnetic units are derivable in two ways : 


via Coulomb's Law for the electric attraction 
(electrostatic units) ; 


via the magnetic analogon of this Law (elec- 
tromagnetic units). 


Accordingly, when both routes are being used, 
a dimensional constant will occur. 
Such a constant is therefore met : 


in Coulomb's Law when using electromagnetic 


units of charge ; 


in the Law for attraction of magnetic poles 


when using electrostatic units of pole 


strength ; 


in Maxwell's Law, which relates the magnetic 
field strength H to the current ¢ producing it, 
when H is expressed in electromagnetic units 
and i in electrostatic units (Gauss’ system of 


units). 


These dimensional constants are equal to the 
velocity of light ¢ or some power of it. They 
could be removed by making the velocity of light 
dimensionally and numerically equal to 1. This 
would mean giving up either the cm or the sec., 
e.g., introducing 3 - 10'° emasa unit of length. 

This might be convenient in the theory of 
relativity, e.g. in describing the energy liberated 
in atomic fission (the energy mc* is now measured 
in units of mass) or in astronomy (light seconds 


simply become seconds), but it is decidedly im- 


practical in general physics. 

A further reduction in fundamenta! quantities 
would be possible by eliminating the gravitation 
constant in Newton's Law for attraction of masses, 
by which, if one sticks to the second, one loses 
the gram too. 

This would leave us with one fundamental! unit 
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and with highly impractical derived units. More- 

over dimensional formulae lose their usefulness, 

since they then provide less insight on relations 

between physical and checks on 

dimensional homogeneity become less powerful 

(note coincidences such as capacity (e.s. egs units) 
length). 

A further difficulty is that dimensional analysis, 
when carried out with too few fundamental units, 
gives too many dimensionless groups. 

A discussion on the latitude in the number of 
fundamental units to be used is given by GUGGEN- 
HEIM [4]. 

Since it is usually avoided to introduce dimen- 


quantities 


sional constants in equations relating quantities 
belonging to the same field of physics, the dimen- 
sional constants, whenever occurring, are normally 
found at the junctions of various fields, i.e., in 
some equations of great fundamental importance, 
such as 


Newton's Gravitational Law | joining kinetics 


Newton's Law of Acceleration | and dynamics 


First Law of Thermodynamics, joining meche 


anics and heat 
Coulomb's Law, joining mechanics and electricity 


Magnetic Attraction Law, joining mechanics 


and magnetism 
Maxwell's Law, joining electricity and magnetism 


The law of gravitation keeps its dimensional 
constant in all systems of dimensions in practical 
use. 

The way of handling the law of acceleration is 
crucial for the dimensions in mechanics and hence 
forms the basis of all dimensions and units in 
general physics and engineering. 

For this reason Newton's Law of acceleration 


will be discussed first. 


Metuops or Usine NeEewrTon’s 
ACCELERATION IN BUILDING a 
OF DIMENSIONS FOR MECHANICS 


4. THe 
Law OF 
SYSTEM 
Newton’s Second Law will be used in the form of 
the statement that force is proportional to mass 
and to acceleration. 


Newton's Law without dimensional constant 


If used without a dimensional constant Newton’s 
Law becomes 


F = ma (static and dynamic systems) (4) 


Since there are obvious stringent reasons always 
to relate the acceleration to the fundamental 
quantities of length and time, there are two 
possibilities left: 


dynamic systems, sometimes called absolute sys- 
tems : 


mass m, length L and time T are fundamental; 
force F is a derived quantity with dimen- 


F = (MLT?).* 


sional formula : 


static systems, sometimes called gravitational, 


technical or engineering systems : 


force F, length L and time T are fundamental; 
mass m is a derived quantity with dimen- 
sional formula m = [FL™ T*). 


Each of these principles can be used with grams, 
kilograms or pounds as mass units or as force 
units, with em, m or ft. as units of length and 
with seconds or hours as units of time, in various 
combinations, 


Newton's Law with a dimensional constant 


If used with a dimensional constant, Newton’s 
Second Law reads : 


g.F = ma (Engl. Engin. system, new version) (5) 


The English Engineering systemt (new version) 
is based on the fundamental quantities of force 
(unit: Ib.,), mass (unit : Ib..), length (unit : ft.) 
and time (unit: either sec. or hr.) with the 
dimensional constant 


g. = [ML T?* F''] (Engl. Engin. system, 
new version) (6) 


Since the same quantity of matter serves to 
define the units of both force and mass, applica- 
tion of eq. (5) to free fall shows that g, is numer- 
ically equal to g. 


* Capital M is used for mass in dimensional formulae. 


+ Name adapted from [1], page 28; however, the 
“ British Engineering system” discussed by CouLson 
and RicHarpson [10] page 2, is a static system. 
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Once the use of four fundamental quantities 
has been decided upon, it is logically impossible 
to relate one of them to the three others. 

The quantities closely associated with force are 
naturally expressed in units derived from Ib. (e.g. 
work, power, pressure, shear stress 7), those more 
closely related to derived from 


mass m units 


ib... (density p). This segregation, however, is 


somewhat arbitrary. Regarding viscosity a 
choice must be made in order that one knows 


dv 
dx 


It is advisable to relate the viscosity 


whether the g. should be placed in + = 9 


or in y pv. 
to the Ib. 


in order also to keep the definition of 
the Reynolds group unchanged. It might be safe 
to distinguish the two definitions by subscripts 
("% and "n° 

If potential energy is expressed in ft.lby and 
kinetic energy in Ib. ft.? sec.*, they will have a 
Also, in the First Law 


of Thermodynamics it should be remembered, 


factor g. between them. 


which of these energy units serves in defining the 
mechanical equivalent of heat. 
These examples can no doubt be multiplied. 
Part II presents an example (number A 10), 


where the unit for pressure is derived from Ib. 


and the unit for shear stress from |b... 

The use of this dimensional system with metric 
units is very rare. E. Scumipr [11] mentions the 
desirability to use it. He uses the symbol kg for 
the kg, kilopond) for 
the kg,. The unit of energy is kpm and not 


only and introduces kp ( 
kg m* sec.*, so that the equality of mechanical 
and electrical units of energy is lost. Since this 
system is not (yet) in use, SCHMIDT carries on in 
technical, i.e. FLT units. 

The English Engineering system (new version) 
has developed from an older version, where the 
two fundamental quantities of force and mass and 
their units were not kept distinguished. 

In this old version both force and mass are 
given the dimensional formula [M] and the unit 
The 
dimensional constant then becomes dimensionally 
and numerically equal to g, i.e. : 


is the pound—without any subscript. 


gF = ma (Eng). Engin. system, old version) (7) 


It is impossible to say that the old version is 


incorrect. It can no doubt be handled correctly. 
However, the use of the same pound for both 
force and mass goes at the expense of clarity, it 
makes dimensional checks and 
analysis less powerful and it must have been a 
source for much confusion. This has no doubt 
induced the development of the new version. 
The old version is still in 


dimensional 


-or recently was 
use by many well-known chemical engineers, e.g. 
Leva [12] and Witne tm [13] in the U.S. and 
Pratt [14] in the U.K. 

The new version is employed at the Mass. Inst. 
Technology (see e.g. McApams [15]), by the 
Dupont Eng. Expt. Station, whose Dr. T. H. 
ardent Shell 
Development Co, 

Waker, Lewts, McApams and GiLuitanp [1], 
page 29, after defining English 
(= absolute), statical ( 
engineering units (i.e. old version) said in 1987 : 


CHILTON is an defender, and by 


dynamical 


gravitational) and 


“ This text will express final equations in engin- 
eering units but wherever ambiguity can arise will 
distinguish between pound mass and pound force. 
It is easier to make derivations in absolute units 
(e.g. mass in pounds, force in poundals) and trans- 
form the final expressions into engineering units.” 


condemnation of the engineering 


It is felt, however, 


Stronger 
system cannot be imagined. 


that by distinguishing lb... and Ib, the greatest 


n 
stumbleblock has since been removed. 

Woke [16] believes to have an argument that 
F = ma is right and g.F 
that the relation differs from all others, which 


only express proportionalities, in that F is ma. 


> ma is wrong. He says 


The tables overleaf gives a survey of some 
important units in the most important systems 
derived according to the four principles described. 


Special Subjects 
The molal units 
In many fields (thermodynamics, kinetic theory, 
stoichiometry) there is a need for a molal mass 
unit. 

It is recommended to introduce [Mol] as an 
additional fundamental quantity in dimensional 
formulae and to regard the molecular weight M 
as a dimensional quantity : 
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M = [M Mol} (8) 


(see footnote § 3 on M in dimensional formulac). 
The units for the quantity Mol are mol, kmol 
and |lbmol. 


Dynamic Systems 





| 
British 
Dynami 
System 


MKS 
gs Practical | 
System System* 


Fundamental 
Units 
length 
mass'!? 
time 
Derived Units 


force 2 


Km cm sec* | kg,,, mi see ba Ib.,, ft. see 
(= dyne) (= newton) | (— poundal) 
| Sen em* sec. ? Kn m* sec.~* Ib. ft.? sec.-* 
(= erg) | (= joule) (=ft.poundal) 
gnem* | kgm a." 
gz, cm *sec.*| kg, m-* see.-* Ib, ft. *sec.4 
dyn. viscosity g,,em ‘see. kg, m 7) Tb. ft.'sec. | 
( poise) | 


work 


density 
sp. weight 





“) ‘The user of a dynamic system will usually omit the 


subscript m. It is added here for extra security. 
{2) Hours may also be used. 


* Warning : The abbreviation “ MKS System ~ 


for both the practical and the technical metric system. 


Static Systems 


Small Metric | 
Static | 





British 
Static 


Technical 
Metric 
Fundamental | 

Units 

length 

force!) 


time 


Derived units | 

gpem™' sec.* | kge m=! sec. | by ft.-' see.” 

(= slug) 
kg,m ft. Ib.» 

kg, m™~* sec.” | Ib., ft. see.” 
kgpm® | Ib ft. 

kg, m™* sec. Ib.¢ ft.~* see. 


Inass 


work | 
density 

sp. weight 
dyn. viscosity gp em sec. 


upom 
£pem * sec.” 


| geem 





is In use 


English Engineering System 





| Old Version 


New Version 
Fundamental Units 
length 
Tass 
time 
force 
Derived Units 
work 
density 
sp. weight 
pressure 
dyn. viscosity, »,,, 
or alternative : n,- ; * sec. 
ily, ft. 


; Ib. gsec. 


dimensional 
constant ‘) 2 





|) The user of a static system will usually omit the sub- 
seript f. 

2) Hours may also be used. 

(3) Note: 


value of g. or g thus becoming 4:17 © 10%, 


Hours are also frequently used, the numerical 


AvoGapro’s number N now becomes a physical 


quantity instead of a pure number : 


6-06 x 10” kmol! 


2-75 x 10% Ib.mol™! 


N 6-06 « 10% mol! 


This procedure is very convenient in dimen- 
checks, For 
Ricuanpson [10], p. 248, give the mean velocity 


' S : , 
of gas molecules as / RT , Obviously a factor 
r 


If R is 
expressed in mechanical work per unit of mass 
(1.2 T*6"'), this 
(ML? T? 


sional instanee, Counsox and 


M has been omitted in the denominator. 


and unit of temperature, R 
omission does not show, but with R 
Mol" @"'} it is found immediately. 
Specific mass (or density) and specific weight (or 
specific gravity) 
If density is denoted by p and specitic gravity 
by y, the hydrostatic pressure at depth A is in 
both the dynamic and the static systems given by 
the two equations 


Pp = pgh (9) 


* R being defined as the universal gas constant. 
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p= yh (10) 


However, the first definition is preferred in 
dynamic systems, where p has a simple meaning 
and the second in static systems. 

Similarly the kinetic energy will be written as 


jor? in dynamic and 4” v* in static systems, 


a 
although these expressions are correct in both 
groups of systems. 

It is that in the 
preferred cquation contains g when gravity has 


thus seen static systems 


no influence and vice versa. 


Specific quantities based on mass or weight 
Attention is drawn to the fact that specific heat, 
entropy, internal energy. etc., are referred to unit 
mass in dynamic systems and to unit foree in 
static systems. 

Thus, these quantities in the two groups of 
differ the 


numerical values are equal and the symbols also 


systems dimensionally, although 
the same. 

Due to this difference in the definitions some 
equations read differently in the two systems. An 
example is the definition of the dimensionless 


group named after Pranpt-. 


nc 
Pr r? (dynamic) 


/ 


Pr ef (static) 


where 
[ML ‘o> 
(QM '@") 
(QL T'e@"'} 


(FL? T) 
(QF"6") 
(QL Te") 
The extra factor is g, not g.. 


5. EXTENSION TO OTHER 


FIELDS OF 


THE 
Puysics 


Heat 


When it is tried to extend the process of deriving 
dimensional formulae of physical quantities to the 
field of heat, the difficulty is met that more than 
one unknown occurs simultaneously so that it is 
impossible to derive the units to be used in this 
field without additional assumptions (ROckER 
[17}). 

In all equations that might serve for such 
derivation heat is equated to a product of tem- 


perature with a quantity such as the gas constant, 
entropy, specific heat or thermal conductivity 
and possibly with purely mechanical quantities 
(mass, area). 

Since there are three unknowns, two assump- 
tions may be made, when the dimensions of the 
third quantity can be found. 

The methods in practical use have onc assump- 
tion in common, Le. temperature is introduced 
as a new fundamental quantity. This is necessary 
in order to be allowed to use the °C or °F. 

The other assumption is : 


Heat is a fundamental quantity. 

Heat is expressed in the units of mechanical 
energy. 

The specific heat is made dimensionless (its 
dimensions are suppressed). 


A suggestion how to derive temperature from 
mechanical quantities was given by RUckex [17]: 
“ The unit of temperature is that necessary to 
produce an increase of Lerg in the mean total 
energy of translation of the molecules contained 
in 1 cm® of gas at 0°C and 1 atm.” Here the two 
assumptions are: heat expressed in mechanical 


units; suppression of the dimensions of R. 


Ad 
The the 


technical metric and various British systems are 


first method is followed when cgs, 
extended with calorie plus °C, keal plus °C and 
BTU (or centigrade heat unit) plus °F or °C 
respectively. The mechanical equivalent of heat J 


is a dimensional constant. 

Ad 2 

The second method is followed when the MKS 
The unit of heat, 


joule, is equal to the mechanical unit of work, 
This method has the advantage 


system is extended with the °C. 


newton metre. 
of not requiring the mechanical equivalent of 
heat. The unit of heat, joule, is also the unit in 
which electrical energy is measured and for this 
reason it is practical to use it in calorimetry. The 
Royal Society has announced [18] preference for 
the joule over the calorie in scientific work. 


Ad 3 


In the third method the specific heat of water is 
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made not only numerically equal to unity — by 
the and BTU — but also 
dimensionally. 


Thus : 


definition of calorie 


heat (M86) 
specific entropy [1] 


(ML! T 


thermal conductivity 


This method is followed by RtCckrer [17], Incer- 
soLt and Zone, [19]. page 6, and by CouLson 
and Ricnarpson [10]. page 3. 

In view of the arbitrariness involved in every 
choice of fundamental! units it cannot be said that 
this method should be incorrect. In fact, making 
the specific heat dimensionless. and unity for 
water at 15°C is not very much more arbitrary 
than making the dielectric constant dimensionless 
and unity for vacuum (in the cgs electrostatic 
units). 

However, it is felt that by this method dimen- 
sional formulae lose much of their clarity to 
workers already used to other systems. While it 
cannot be denied that the method may be useful 
in heat transfer, it certainly is in the way if we 
see the use of one coherent system of units for 
the entire field of physics and engineering as the 


ultimate aim. 


Electricity and magnetism 

Here, too, the difficulty arises that the units 
cannot be derived from purely mechanical units 
without an additional assumption, Coulomb's 
Law, for instance, containing two electrical 
quantities, charge and dielectric constant. 


The two ways followed are : 
suppression of the dimensions of one electrical 
quantity : 
in cgs e.s. system (dielectric constant), 
in egs ¢.m. system (magnetic permeability). 
making one electrical quantity fundamental : 
in practical system. 
For the cgs system reference is made to the 
discussion, p. 134. 
In the practical system both the dielectric 


constant « and the magnetic permeability » are 
no longer pure numbers. Each electrical unit can 


be derived in one manner only. No electrical unit 
can be derived from mechanical units only (see 
Pout [8)}). 

Coulomb’s Laws for electrostatic and for mag- 
netic repulsion in vacuo in the practical system 
are written 


®, o, 
+ mg" 


e; C2 


+ we,’ 


with ¢,, ¢g, electric charges, and ¢,, ¢,, magnetic 
fluxes. 

Whether «, and yy are interpreted as dimen- 
sional constants or properties of the vacuum is 
probably a matter of taste. 


€) = permittivity = dielectric constant of empty 


space 


107 A sec. 9-888 x 10" A sec. 
tr (3 x 10°)? Vm welia Vm 


Hy = permeability = induction constant of empty 
space 
, V sec. 
tr 10°° 
Am 


with 
1 


£o Ko 


(3 x 10° m/sec.)* = c* 


The factor 47 is introduced for the purpose of 
rationalization (see page 133). 

From the user’s point of view it is immaterial 
which extra electrical quantity is used in fixing 


the practical units. Proposals are : ohm, ampére 


~ V sec. 


(actually preferred) and yp, = 47 x 10° 


Am 

Dimensional formulae and units, expressed in 
fundamental units may cumber- 
some. WILLIAMS [20] in a treatise on the practical 
system, e.g. states the unit of «, to be m™* kg™ 
sec.* A? 


become very 


Units, the physical significance of which is easy 
to grasp, are obtained when expressing them in 
terms of V, A, m and sec. [21, 22]. 


The units for ¢, and yy as given above are seen 
to be symmetrical in structure and convert into 
each other if V and A are interchanged. This very 
convenient symmetry is found throughout electro- 
statics and magnetism. 
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Compare : Unit Unit 
V — V sec. 
electric field strength E magnetic induction, B ; 
m m 
: A sec. , V sec. 
dielectric displacement, D ; inductance, L 
m A 
capacity, C A see. For further literature on the practical system 
with ' \ reference is made to: Van Dick [22], Nether- 
A lands Standardization Sheets [21], Pons. [23], and 


agnetic field strength, H ; 
catinaieitien . m Corneuivs and Hamaker [5]. 
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On the steady state fractionation of multicomponent and complex 
mixtures in an ideal cascade 


Part 3. Discussion of the numerical method of calculation 
AnprReEAS Acrivos* and Neat R. AMuNDSON** 
( Received 26 January 1955) 


Summary—In this paper the formulae necessary for rectification calculations on complex 
mixtures are summarized and developed for numerical computation. These formulae involve 
the computation of improper integrals for which a method of calculation is developed. A numerical 
exumple is considered in some detail since these calculations are not of a routine nature. The 
calculations were carried out in the I.B.M, 602 A. Calculating Punch. The plate by plate method 
is also considered. 


Résumé — ans cet article, les auteurs ont résumés et développés, en vue d’applications numériques, 
les fornules nécessaires au calcul de la rectification dans le cas de mélanges complexes. Ces for- 
mules impliquent le caleul d'intégrales impropres pour lesquelles ils exposent une méthode de 
calcul. Ils étudient en détail un exemple numérique puisque de tels calculs ne se rencontrent pas 
couramment. Les calculs ont été publiés dans le * 1.B.M. 602 A. Calculating Punch.’ Les auteurs 
considérent aussi la méthode * plateau par plateau.” 


RESUME AND FurTHER DEVELOPMENT A mass balance around the top n stages gives 


oF NECESSARY FORMULAE : 

. t Lol p) , 

In the previous two parts the general theory of Yn +i(P) R4 j ta P) "RoW (1) 

the distillation of complex mixtures has been 

developed for an ideal cascade. Formulae were while the equilibrium relation is assumed to be, 
developed, under the usual simplifying assump- for continuous mixtures, 

tions, which enabled one to make computations Punt P) 
@ 


involving the number of stages required or the <,\P) 


compositions attained by either a plate by plate ' py, p) dp 


method or by means of a formula in closed form 
: ri 


involving some improper integrals. In order that 
In Eq. (1), 2o(p) represents the distribution 


reference to the previous papers be minimised the 
function of the overhead product. Similarly, for 


necessary formulae are summarized here. 
i ; the stripping section, 
If z,(p) and y,(p) are the distribution functions a it ila 
of the liquid and vapour, respectively, from the *(p) 
) = 
nth stage, then Yo R* + 


| =e) dp = 1 z,*(p) = 
i | py,*(p) dp 
| yp) dp = 1 . 


0 


PYys*(P) 
_£ 


where now y*,(p) and 2,*(p) are, respectively, the 


* Department of Chemistry and Chemical Engineering, University of California, Berkeley, California. 
** On leave at Department of Chemical Engineering, Cambridge University, England. 
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distribution functions of the vapour stream and 
R* is the 
reflux ratio in the stripping section, and x9*(p) 1s 
the distribution function of the bottom product. 


It has been assumed also, for simplicity, that the 


the liquid stream leaving the sth plate, 


relative volatility of every component is the same 
in both sections of the column, but this is not at 
all necessary. 

Equations (1) and (2) are the basi 
and in Part I it was shown that the ec 


equations 
ym position 


on any plate is given by 


UA Pp) 


symbol P in front of 


integral denotes its Cauchy principal value, 


improper 
that 


where the 


P E To\ q) d 


r 
A—q / 


. 
“o 


hq. (5) then expresses the functional dependence 
of x,(p) on both p and n in a closed form, It 
enables one to calculate vr, P) directly, that is 
without making use of any intermediate steps, if 
the distribution function of the overhead pro- 
duct, zo(p) is known. Moreover, it can easily be 
shown [1] that Eq. (5) reduces to FENSKE’s 
well-known equation as R + @ 

Eq. (5) can, however, be put into a more con- 
venient form by introducing a new function, 
G(A), defined as, 

C(A) (A\*? 
(;) 


Bia) (9) 


G(A) = 


where n = N denotes the feed plate and h is 


some point on the p-axis which is, as yet, arbi- 
trary. Then, at the feed plate 


p{ 0), 
ptlp) % 
R @ 
G(A) da 


G 
1A 4 \P) w(p) 
Pro P) 


* NP) 


vu 
\ useful asymptotic expression can be obtained 
from the above, since, when 2,(p), for a given 
value of p, is small enough that 


x 


" GA 
p tol P) P \ \ da 


Pp 


G(p) w(p) 


then 
_ Gp) w(p) 


R|G(A) da 
Finally, the point h should be chosen in such a 


as to facilitate 
In other words, h should be « hosen 


way the numerical evaluation of 
the formulae. 
so that 

G(p) « for p 
and 


Gp) 10 for p> h 


Entirely analogous formulae are obtained for 


the stripping section. It is found that 


r,*(p) 


Pp 7 'C*(A) 

prop) 3 2” 
k* = 
= C*A) dA 


_?P 'C*(p) 


1A w*(p) 


0 
where C*(A) and w*(A) are defined by expressions 
analogous to equations (7) and (8), respectively, 
that is 
A r*(A) 


C*(A) = [w*(A)]}* j [mA q*(A)}? 


(14) 


(15) 


~ * 

= R* — P 9*o ) ag 
. q — A 

0 
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Again, in order to facilitate the numerical work, 
define 


C#(A e+e 
Ga) =~! (| (16) 


~ C*(2) 


where the symbol s S denotes the feed tray 
and / is somewhat arbitrary but must be chosen 


Then 


properly. 


ry*(p) 


P [ ot) ain 
pry*(p) 3? 
R* ; 
(i*( A) dA 


G*| 
ete w*( p) 
w..*(p 

Pla \P) 


(17) 


vu 


small cnough 


and if a,*(p) is 


G*(p) w*(p) 
rs*(p) -> ——@ (18) 





R* | G*(A) da 


0 


Moreover, at the feed plate 


rg*(p) ry(p) 
and therefore 
GF Gr 
P| ar ' ae w*(p) 
Tol P) ZR 0 “4 A pP ql oA P) (19) 
r . ?) k* ©, : k . 
0 (7 P| G(A) G(p) . ) 
. A—p PrP) 
where 
| G(A) dx 
Z=°— (20) 
| G*(A) da 
0 


If either z,*(p) or a(p) is small, then Eq. (19) 
can be simplified to give, respectively, 


G*(p) = 


R* p2o{p) f. G(A) G(p) 
“Eo Pi —— di + 21 
ZR w*(p) A—p P XP) seid 


ri 
and 
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G(p) 
RZ pro"(P)| p (GA) gy , EXP) yy] (29) 
R*u(p) Jp—aA p Xy*(p) 
Tuk NuMmMeErRIcCAL EvAaLuaTION 


OF AN IMPROPER INTEGRAL 


Since all the important formulae of the preceding 
section included the Cauchy principal value of an 
improper integral, it was thought that a method 
should be pre sented here for numerically calcu- 
lating such a function. Let us suppose then that 


f(r) is a single-valued, integrable, continuous 


function of x. which vanishes everywhere for 


Then 


negative a 





J ) y . ; ) \ 
I(A) / hd KE ; A dx (28) 
A—J , I 
or 
IA) (A, x) da (24) 
*o 


The above integral can now be evaluated numer- 
ically in a number of ways, but use will be made 
here of Simpson’s rule, that is, if the x-axis is 
subdivided into m equal intervals of length A, 
where m is even, then according to Stmmpson’s 


formula 


° m/2 
Kay =™ je ) 42S $(A, 2kh — h) 
3 2 1 
m/2-1 
+ & ¢$(A, 2kh) 4 aati] (25) 
l ~ 


The various ordinates can, of course, be easily 


evaluated numerically except for ¢(A, 0). How- 
ever, 
$(2,0) = lim LA—2) —SAF2) _ _g ry) (96) 


r>0 a 
A convenient method for obtaining the derivative 


of a function has been described by RuTLEDGE [2], 
and below are two of his results : 


(a) If, for 2 — 2h Sq Sa + 2h, f(g) can be 
represented as a polynomial! of the fourth degree, 
then 
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1 
I) = tay [see — 2h) - Fe + 2h) + 


+ 8{f(@ +h) —f(@—h)}} (27) 

(b) If, for z — 3h Sq Sa + Bh, f(q) can be 
represented as a polynomial of the sixth degree, 
then 


: ie ; 
St’ (@) ia0h 3 [f(# 3h) — f(a — 8h)] 
+ 27 [f(x — 2h) — f(x + 2h)] 
+ 185 [f(a +h) —f(a — h)) | 


Finally, define 


f(A — kh) — f(A + kh) 


z (A, kh) (29) 


Table 1. 
PprpAp) = 2-25676 p(p 


Numerical values for 


le 





pry) 


PrrP) 


prep) 


0-04153 


0-03325 


56215 
490554 
41999 
33798 
25112 
16120 
O60B8S 
‘OTS67 
SSS04 
80182 
71830 
63913 
56490 
49601 
43268 
37502 
32327 
0-27637 
0°23503 
0-19863 
0-16683 
0-13927 
0-11556 
0-09530 
0-07813 
0-06367 
0-05158 


0-00000 
0-00591 
0-02458 
0-05709 
0-10408 
0-16563 
0-24132 
0-33019 
43077 
54117 
65909 
78193 
90691 
03112 
15167 
26579 
37086 
46460 
54507 
61073 
66043 
69361 
71001 
“70989 
69391 
1 -66306 
1-61862 


0-02048 


to 
= 
oo 


& 


0-02095 
0-01649 
O-OLZFG 
0-01003 
0-00776 
0-00596 
0-00456 
0-00347 
0-00261 
0-00197 
0-00147 
0-00110 
0-0008 1 
0-00060 
0-00043 
0-00032 
0-00023 
0-00016 
0-00012 
0-00008 
0-00006 
0-00004 
0-00003 
0-00002 


ot 0 tO OW OW WS 8S 


eee eee 








Then 


2 m/2 
(A) = 5 [As +2 2D 2(A, 2kh — h) + 


m/2-1 
t+ 2 
1 
The above is believed to be the simplest form 
into which Eq. (23) can be reduced. 


1 
z(A, 2kh) + 


z (A, mh) 
- (30) 


A NUMERICAL EXAMPLE 


Let z,{p), the distribution function of the feed, 

be given by 

. 4 2 (p-1)* 

x Ap) (p—1pPe forallp>1 (81) 
Vr 

The vapour to liquid 

Moreover, it is 


and 2,(p) = 0 for p < 1. 
ratio in the feed is equal to 1. 
specified that 


2,(1-50)/x9*(1-50) = 20-0 


x*,(2-95) ‘x9(2-95) = 20-0 (82) 


and that R 2-15. Tables 1 and 2 show some 


numerical results. 
It is found from the calculations that 


D = 0.509138, RD = 1.09463 and R*W = 1.59462 


Table 2. A preliminary estimate of the functions 


A Dad) and AWz,*(A). 





DaglA) | AWag*(A) 


| ADaglA) | AWag*(A) 
0-603892 | 1-01470 
| 053824 | 102391 
047316 | 1-02218 
0-41089 | 1-00011 
0-35162  0-98636 
0-20644 | 0-95468 
024562 | 0-91558 
| 0- 20022 0-86965 
015049 | 0-81918 
| 012386  0-76508 
0-09312 | 0-70870 
0-06714 | 0-65116 
0-04495 | 0-50419 
0-02753 | 0-58737 
001550 | 049200 
 0-00420 | 0-42848 


0-00350 
0-03566 
0-07803 
0-12919 
0-18832 
0-25574 
0-32971 
0-40768 
0-48836 
0-57058 
0-65080 
0-72739 
0-79864 
| 0-86225 
| 0-91650 
| 0-96073 
0-99355 


0-537 67 
0-62343 
0-70390 
O-77772 
0-84281 
0-89593 
0-93608 
0-96318 
0-97624 
0-97454 
0-95993 
| 0-93304 
0-89497 
0-84776 
| 0-79339 
0-73318 
0-66952 
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Table 3 (a) 





A 
45625 
4625 
46875 
“4750 
“48125 
“4875 
49375 
500 
“50625 
51250 
“51875 
“5250 
‘53125 
53750 
54375 
5500 


eet 


w*(A) 
0-01359 
0-01078 
0-03448 
0-05752 
0-07989 
0-10160 
0-12264 
0-14301 
0-16100 
0-17881 
0-19645 
0°21391 
0-23120 
0-24832 
0°26526 
0- 28202 


AW x9*(A) 


0-00696 
0-01058 
0-01436 
0-01830 
0-02240 
0-02666 
0-03108 
0-03566 
0-04048 
04545 

0-05052 
0-05575 
0-06111 
0-06661 
0-07225 
0 07803 


c*%A) 
1038-61296 
85-50896 
43°-95479 
27°30912 
19-50499 
1517744 
12-48312 
10-66428 

“49112 
56595 
‘81731 
‘19858 

+ 67832 

}° 23433 





60625 
6125 
“61875 
6250 
63125 
“63750 
64375 
“6500 


we ee ee ee ee ee ee en ee 


mw*( A) 


0-29922 
0-31607 
0-33258 
0-34873 
36453 
37998 
“39508 
40983 
42413 
43811 
45177 
046510 
047830 
0:49096 
050827 


0-5152 


AWa*o( A) | 
0-08400 
0-09007 
0-09628 
0-10262 
0-10907 
0-11565 
0-12236 
0-12919 

13613 

14319 
15039 
‘15772 
16517 
‘17276 
18047 
18832 


C*(A) 
5°20845 
4°93945 
470195 
4:49048 
4°30086 
4°12974 
3-97437 
3°83259 
3-70342 
3-58413 
3-47351 
3°37055 
3-27328 
3°18346 
3-09915 
3-01980 





Table 3 (b) 





2-8125 
2-825 
2-8375 
2-850 
2-8625 
2-875 
2-8875 
2-00 
2-9125 
2-925 


wi A) 


0-34350 
0-32334 
0-30312 
0- 28284 
0-26256 
0°24218 
022169 
0-20111 
0-17978 
0-15878 


ADwof A) 


0-08627 
0-07 966 
0-07328 
0-06714 
0-06115 
0-05545 
0-05005 
0-04495 
0-04032 
0-03588 


CA) 


444756 
4-70278 
499104 
532286 
5-70190 
6°14927 
668687 
734649 
8°22783 
0-33889 





30125 
8-025 
83-0375 


wi A) 


-~QO- 13811 
0-11777 
0-099040 
0-08027 
0-06037 
0-03971 
0-01637 
0-00646 
0-02878 


ADaof A) 


C(A) 





0-03161 
0-02753 
0-02445 
0-02142 
0-01844 
0-01550 
0-01261 
0-00976 
0-00696 


10°78167 
12-72591 
15°29186 
19-26582 
25-99771 
38-75103 
67 -74756 
08-11714 
52-58279 





Table 4 (a). 


Numerical values for G(A). 


= 86 and 





2-70 
2-7125 
2-725 
2-7375 
2-750 
2-7625 
2-775 
2-7875 
2-80 
2-8125 
2-825 


G(A) 


GA) 


G(A) 


G(A) 





0-01082 
0-01476 
0-01820 
0-02119 
0-02380 
0-03273 
0-04057 
0-04743 
0-05338 
0-06575 
0-08120 





2-8375 
2-850 
2-8625 
2-875 
2-8875 
2-90 
2-9125 
2-925 
2-9375 
2-950 
2-9625 


0-10059 
0-12510 
0-15619 
0-19619 
0-24833 
0-31735 
0-41317 
0-51768 
0-73023 
1-00000 
1-41732 





2-975 
2-9875 
3-00 
83-0125 
8-025 
3-0375 
3-050 
3-0625 
3-075 
3-0875 
3-100 


2°10394 

334344 

586495 
12-05879 
20°52533 
12-91910 
3-47029 
2°29453 
1-70332 
134813 
111177 


31125 
3-125 
3°1375 
3°150 
3-1625 
3175 
3-1875 
3°20 
3-30 
3-40 





0-95156 
0-82655 
0-72609 
0-64354 
0-57682 
0-51931 
0-46889 
0-42494 
0-21647 
0-12006 
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Table 4 (b). Numerical values for G*(A). N 


36 and S 





G*(A) 
200 0-16677 
‘250 0-30875 
300 0-54662 
‘3125 0-63335 
325 0-73229 
3375 0-84708 
350 0-98126 
‘3625 115219 
375 1-36671 
‘3875 1-64204 
Seal 201011 
4125 2-60980 


GA) 


364364 
6-00723 
13-76852 
49-49106 
32-27127 
13°12046 
645393 
3-65310 
2°25500 
1-47272 
1-0000 
O-TOBO2 


— eee eee et 








G*(A) GA 


0-01713 
0-01328 
0-010383 
0-00805 
0-00629 
0-00493 
0 -00387 
0-00304 
0-00240 
0-00189 
0-00150 


0-30895 
0-37006 
0-27198 
0-20145 
0-15033 
0-11295 
0-08519 
0-06449 
0-04905 
0-038749 
0-O02877 
0-02216 


63125 
6375 
64375 
“6500 








Next, choose A 2-95 and / 1-50 and calculate 
the functions w(A), C(A), GA), w*(A), C*(A) and 
G*(A). It was found by repeated trial and error 
that letting N = 36 and S§ 


factory. Many of the values shown in Tables 3 (a) 


54 was fairly satis- 


and 3(b) were obtained by a parabolic inter- 
The values of G(A) and G*(A), 
in Tables 4(a) and 4(b) were either calculated 
from Eq. (9), 2-70 = A s 83-0875, and Eq. (16), 
1.45624 A 1.6500, respectively, or from Eq. 
(22), A> 3-05, and Eq. (21), A 1-450, respec- 
The tabulated values of G(A) and G*(A) 
for the latter two ranges, respectively, are only 
One from 


polation. shown 


tively. 


approximately correct. also finds 


Eq. (20) that 
Z 0.86938 


Finally, by means of Eq. (19) and the obvious 
relation 


rp p) = Dap) Wr,*(p) 


one calculates numerically the functions AD a (A) 
and AW 2,*(A). 
Table 5. 

A comparison now between Tables 5 and 2 will 
reveal that the two functions ADz,(A), that is, the 


The results are presented in 


one which was assumed, and the one which was 
calculated, are not very different. The agreement 
between said two functions however cannot be 
considered, in this particular problem at least, 
satisfactory, because, as can be verified by a 
numerical calculation, a small change in w(A) and 


Table 5. Calculated numerical values for ADa¢{ A) 


and AW2,*(A) 





ADagfA) | AWa*of A) AWay*(A) 


ADaof A) 


1-04581 
1-05433 
1-05109 
1-03700 
1-01289 
0-97981 
0-93918 
0-89229 
0-83841 
O-78113 
0-72133 
0-6603 1 
0-59871 
0-53737 
0-48021 
0-42848 


0-53767 
0-62343 
0-71266 
0-787 52 
0-86159 
0-88949 
0-92525 
0-947 36 
0-95564 
0-95058 
0-93280 
0-90327 
086390 
0-81578 
0-76080 
0-70081 
0-63765 


0-003850 
0-03556 
0-6927 
0-11939 
0-17054 
0-26219 
B4054 
$2350 
50896 
SO4AY 
0-67793 
O-75716 
082971 
0-89423 
0-94909 
0-99310 
1-02541 


38290 
32509 
27131 
22202 
‘17759 
14026 
“10781 
0-08049 
0-05799 
0-04042 
0-02753 
0-01580 
0-00420 


ws w Ww 


te 








w*(A) will alter the functions C(A) and C*(A) con- 
siderably, especially in the sensitive regions 
2-95 « < 3-10 and 1-40 <A< 1-50. This is 
due to the fact that R was deliberately chosen 
close to Rains 
rate, the numerical values of both N and S would 
markedly depend on the assumed distribution 
functions of the two products. This of course, 
would be expected; on the other hand, it is 
somewhat surprising to find that the functions 
C(A) and C*(A), and therefore G(A), G*(A), N and 


and therefore, in this case at any 
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S, do depend in a very pronounced manner on 
the assumed values of the distribution function 
for 2:95 < A < 8-05 and 1-45 < A < 1.50, respec- 
tively. It is true that ADa, (A) is small for 
2:95 < A < 8-05, yet it is precisely the value of 
ADz,{A) in the above mentioned interval which, 
together with w(A), primarily determines the 
function G(A); entirely analogous considerations 
can be made concerning AWz,*(A) for 1.45 < A < 
1-50. We see then that G(A), which is large only 
for 3-05 > A > 2-95, depends markedly on the 
whole distribution function Dz,(A), and not only 
on part of it, as one would have expected. So, 
one of the difficulties which one encounters in such 
problems is to be able to estimate the distribution 
functions of the products in every detail. That is 
even in the regions where they come numerically 
very small. And it is because this estimate of 
the distribution 
rapidly that the analytic form of the solution is 
preferable to that obtained by the iterative 
Moreover, by means of the analytic 


functions can be made more 


method. 
solution, one can see more clearly, than by the 
iterative solution, which factors in any numerical 
problem primarily affect the answer and to what 
an extent. 

This numerical example, then was presented 
in order to illustrate the application of the formu- 
lae to an actual problem, to point out certain of 
the difficulties encountered in such problems, and 
to bring out certain of the characteristic features 
of the analytic solution which should make it 
superior to that derived by the iterative method. 
In order to solve this problem completely one 
should revise the original estimate on the function 
ADa,{A) and repeat the procedure outlined in this 
chapter, until the assumed and the calculated 
values for ADz,(A) become essentially identical ; 
this, however, shall not be attempted here. 
PLaTteE METHOD 


Tue PLATE BY 


In Part I it was shown that the composition of 
the liquid on any tray is given by 


| a, X,, dp 
0 
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where 

X(p) = RpX,(p) + p | 2,(p) X,(p) dp 
n+l - 

with 

Pp 


X,(p) R 


Successive substitutions will show that if M, is 
defined by 


M, To( p) p* dp 
then 
X,{ p) p* pt 4 M, 
‘ ; ‘ | ‘ 
X,(p) = Rp® + p*M, +p | M, + ,M 
and, in general 
X,(p) = P,(p) 


where P,(p) is a polynomial of the nth degree 
having only positive coefficients. Therefore P,(p) 
is positive and strictly monotone increasing. 


Further 


with 


*8 


N rol p) P,(p) dp 


* 


0 


X,, has a form independent of the distribution 
function and the same may be said for N,. Also 
the coeflicients in P, depend only upon the reflux 
ratio and the successive statistical moments of 
the distribution function. The method of calcu- 
lating the composition on a particular plate, 
knowing the composition of the overhead is then 
a straightforward task on a modern computing 
machine. A programme for doing this has been 
prepared for the 1.B.M. 602 A. 


CONCLUSIONS 


The solution just outlined to the basic equation 
for the rectification of complex mixtures in an 
ideal plate column has all the well-known 
attractive features of an expression in a closed 


form. At this point, however, one might be 
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tempted to ask whether or not it would have been 
simpler to break up the complex mixture into a 
finite number of distinguishable components and 
then make use of the equations already developed 
by Unperwoop and the other investigators, 
rather than employ the formulae of this paper. 
On the contrary, it would be preferable, in 
general, to replace a large number of distin- 
guishable components b ya continuum, because 
the expressions developed here are easier to 
handle than the corresponding ones for the finite 
problem. For, it will be remembered that in 
UnpDERWoOop’s method one has to calculate by 
trial and error some of the roots of the equation 


E PU)*0) _ p 
jar A— PJ) 


where m is the total number of components in 
the mixture, whereas in the method just outlined 
no such trial and error calculation is necessary. 
And at any rate, the expressions presented in this 


paper can be easily evaluated numerically by 
means of a digital computer, a feature which 
should add, naturally, to the usefulness of this 
analytic solution. 


NOTATION 


z,(i) = mole fraction of the ith component in the liquid 
stream leaving the nth plate in the column 

y,(t) = mole fraction of the ith component in the 
vapour stream leaving the mth plate in the 
column 

1/a = the reciprocal of the relative volatility 

«,(p) = the distribution function of the liquid stream 

leaving the nth plate 

the distribution function of the vapour stream 

leaving the nth plate 


4,(P) 


to(p) = the distribution function of the overhead pro- 
duct 
C(A), wt A), and G(A) are functions defined, respectively, by 


a Eqs. (7), (8) and (9). 


p * f(z) 


—s dz is the Cauchy principal value of the improper 


‘ integral 
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